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probability  (MAP)  is  derived  for  image  signals  detected  at 
low  light  levels.  These  signals  suffer  from  Poisson  noise 
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and  blurring  degradations. 

The  low  level  photon  resolved  image  signal  is  modeled 
as  an  inhomogeneous  Poisson  point  process.  The  photon 
noise  is  inherent  in  any  detected  image,  and  is 
particularly  serious  at  low  light  levels.  At  these  low 
light  levels,  the  emission  of  photons  is  described  by  a 
Poisson  point  process,  with  the  average  rate  of  emission 
proportional  to  the  integrated  intensity.  The  blurring 
degradation  model  in  the  system  includes  space-variant  and 
space- invar iant  effects  such  as  atmospheric  turbulence, 
linear  motion,  diffraction,  and  aberrations.  The 

estimation  is  performed  assuming  that  the  photon  events 
i counted  in  each  detector  are  independent,  Poisson 

distributed  random  processes  for  the  large  time-bandwidth 
product  case.  Since  the  variance  of  the  Poisson 

^ distribution  is  identical  to  its  mean,  the  Poisson  noise  is 

neither  multiplicative  noise  nor  a linear  additive  Gaussian 
noise,  and  is  generally  signal-dependent.  It  has  been 
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model  has  improved  performance  because  the  MAP  filter  can 
be  generalized  to  linear  or  nonlinear  image  models  and  to 
noise  models  different  from  additive  Gaussian  noise.  In 
addition,  the  MAP  filter  can  be  a local  adaptive  processing 
filter  and  extended  to  the  case  of  space-variant  blurring. 
It  also  has  been  shown  that  image  models  with  a 
nonstationary  mean  and  stationary  variance  give  useful  a 
pr ior i information  for  the  MAP  filter.  The  MAP  estimation 
equations  are  nonlinear  and  have  large  dimensionality.  A 
sectioning  method  with  a Newton-Raphson  solution  has  been 
adapted  to  cope  with  these  problems.  It  has  been  shown 
that  the  strategy  is  an  effective  and  fast  way  to  solve 
nonlinear  MAP  estimation  equations . 


The  Cramer-Rao  lower  bound  (CRLB)  on  the  mean-square 
estimation  error  of  the  MAP  unbiased  estimate  is  derived 
for  the  Poisson  noise  model.  It  is  shown  to  be  a very 
useful  bound  for  finding  the  best  suboptimal  sectioning 
filters. 


Finally,  a comparison  between  the  performance  of  the 
MAP  filter  and  that  of  the  linear  minimum  mean-square  error 
( LMMSE ) filter  is  made  for  Poisson  noise  models.  The 
performance  of  the  MAP  filter  is  much  better  than  that  of 
the  LMMSE  filter.  The  LMMSE  filter  works  very  well  for 
higher  signal-to-noise  ratios,  but  the  MAP  filter  works 
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CHAPTER  1 


’I 
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. 4 


INTRODUCTION 

.1.1  Introduction 

Image  restoration  can  be  viewed  as  an  estimation 
process  in  which  operations  are  performed  on  observed  or 
measured  noisy  data  to  estimate  the  object.  More  clearly, 
image  restoration  is  the  estimation  of  the  original  image 
signal  by  both  blur  removal  and  noise  suppression.  Image 
enhancement  is  the  attempt  to  improve  the  appearance  of  the 
image  for  better  human  viewing  or  machine  processing. 
Hence,  image  enhancement  may  not  specifically  need 
knowledge  of  the  degrading  phenomena.  However,  in  order  to 
effectively  develop  an  optimal  restoration  filter  with 
various  criteria,  it  is  necessary  to  characterize 
quantitatively  the  image  degradation  effects  of  the 
physical  image  system.  Image  restoration  begins  with  a 
model  of  degradation  effects,  assumes  given  a pr ior i 
information  and  then  develops  an  optimal  filter  to  obtain  a 
restored  image.  Hence,  accurate  image  modeling  and  more  a 
pr ior i information  are  often  the  key  to  effective  image 
restoration. 
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■ifn  this  dissertation,  optimal  restoration  filters  are 
developed  in  the  sense  of  maximum  a posterior i probability 
(MAP)  and  maximum  likelihood  (ML)  for  blurred  image  signals 
detected  at  low  light  levels.  This  low  level  photon 
resolved  image  signal  is  modeled  as  an  inhomogeneous 
Poisson  point  process.  The  photon  noise  (which  we  call 
Poisson  noise  throughout  this  thesis)  is  inherent  in  any 
detected  image  signal  particularly  at  low  light  levels.  At 
these  low  light  levels,  the  emission  of  photons  is 
described  by  a Poisson  process  with  the  average  rate  of 
emission  proportional  to  the  integrated  incident  intensity. 
The  estimation  is  performed  assming  that  the  number  of 
photon  events  counted  by  the  detectors  are  independent, 
Poisson  distributed  random  processes  for  a given  unknown 
object  radiance.  Since  the  variance  of  the  Poisson 
distribution  is  identical  to  its  mean,  the  Poisson  noise  is 
neither  multiplicative  noise  nor  linear  additive  Gaussian 
noise.  It  is,  indeed  a signal-dependent  noise. 


1 . 2 Organization  and  Contributions  of  the  Dissertation 


In  next  chapter,  we  discuss  three 
some  system  models  for  image  noise 
inhomogeneous  Poisson  process  model 
counting  system.  In  Chapter  3 , we  r 
linear  and  nonlinear  image  restoration 
noise  models  and  their  motivation 


image  models  and 
; we  also  present  an 
which  is  a photon 
eview  some  important 
filters  for  Poisson 
for  the  work  in  this 


dissertation.  We  also  review  other  restoration  filters  for 
different  image  noise  models.  In  Chapter  4,  we  develop  and 
implement  an  MAP  filter  without  blurring  degradations  for 
the  Poisson  noise  model.  In  Chapter  5,  we  develop  and 
implement  MAP  filters  with  blurring  degradations  for  the 
Poisson  noise  model.  In  Chapter  6,  we  derive  Cramer-Rao 
lower  bounds  (CRLB)  on  the  estimation  error  for  MAP  filters 
and  discuss  the  results.  In  Chapter  7,  we  compare  the 
restored  image  performance  of  the  MAP  filter  with  that  of 
the  LMMSE  filter.  In  Chapter  8,  we  conclude  this 
dissertation  and  discuss  future  research  on  the  problem. 

The  specific  research  contributions  of  this 
dissertation  are  now  summarized.  A model  for  photon 
resolved  low  light  level  image  signals  detected  by  a 
counting  array  is  developed.  These  signals  are  impaired  by 
signal  dependent  Poisson  noise  and  linear  blurring. 

An  optimal  restoration  filter  based  on  maximizing  the 
a posterior i probability  density  (MAP)  is  developed.  A 
suboptimal  overlap-save  sectioning  method  using  a 
Newton-Raphson  iterative  procedure  is  used  for  the  solution 
of  the  high  dimensionality,  the  nonlinear  estimation 
equations  for  any  type  of  space-variant  and  invariant 
linear  blur.  An  accurate  image  model  with  a nonstationary 
mean  and  stationary  variance  is  used  to  provide  a pr ior i 
information  for  the  MAP  restoration  filter.  The  Cramer-Rao 
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lower  bound  (CRLB)  of  the  unbiased  MAP  restoration  filter 


U 

• j 

is  derived.  Finally,  a comparison  between  the  MAP  filter 
and  a linear  space- invar iant  minimum  mean-square  error 
(LMMSE)  filter  is  made. 


CHAPTER  2 


SYSTEM  MODELS  FOR  IMAGE  NOISE  AND  IMAGE  MODELS 

2 . 1 Introduction 

Using  a statistical  approach  to  develop  an  optimal 
restoration  filter,  the  effectiveness  of  the  algorithm 
depends  on  the  completeness  of  the  statistical  description 
of  image  and  noise.  Hence,  we  investigate  image  and  noise 
models  in  detail  in  the  following  two  sections. 

In  section  2.4,  we  present  an  inhomogeneous  Poisson 
process  model  (we  refer  this  inhomogeneous  Poisson  process 
as  the  Poisson  noise  degradation  throughout  this  thesis) 
which  is  a photon  counting  systtem.  In  section  2.4.1  we 
investigate  a general  photon  counting  system  containing 
blurring  and  Poisson  noise  degradations.  In  section  2.4.2, 
we  first  discuss  a single  detector  model  and  derive  its 
statistical  properties.  In  section  2.4.3,  we  extend  the 
single  counter  to  an  array.  In  section  2.4.4,  we  discuss 
quantum  limitations  of  photon  resolved  image  signals,  and 
in  section  2.4.5,  we  discuss  the  simulation  of  images  with 
Poisson  noise  at  different  ensemble  mean  signal-to-noise 
ratios.  A comparison  is  made  between  Poisson  noise 
degraded  images  and  linear  additive  Gaussian  noise  degraded 


images.  Finally  some  conclusions  of  this  chapter  are 
presented . 


2.2  System  Models  for  Image  Noise 

2.2.1  Linear  Additive  Gaussian  Noise  Model 


This  model  is  most  often  used  for  image  formation  in 
the  field  of  digital  image  processing.  Its  block  diagram 
is  illustratea  in  Fig.  2.1,  where  the  image  g(x,y)  and 
object  f(x,y)  may  be  considered  intensity  functions  of  two 
spatial  dimensions  (x,y),  and  h(x,y)  is  the  point-spread 
function  (PSF)  or  impulse  response  of  the  imaging  system. 
Because  the  linear  blurring  degradation  in  all  image  noise 
models  is  the  same  as  in  this  model,  we  investigate  more 
details  about  the  PSF  in  this  section.  The  function  n(x,y) 
is  additive  noise  which  is  signal-independent  and  Gaussian. 
To  unify  the  notation,  we  denote  continuous  functions  with 
(x,y)  and  discrete  functions  with  (i,j)  throughout  this 
thesis.  The  mathematical  representation  of  Fig.  2.1  is 


g (*,y) 


oo 

h(e,n;x,y)f(e , n) dedn+n (x,y ) . 

— oo 


(2.1) 


This  equation  is  a first  order  Fredholm  integral  equation 

plus  a random  noise  component,  where  (e,n)  is  the  spatial 

coordinate  of  the  object  of  interest  and  h(e,n;x,y)  is  a 

general  space-variant  point-spread  function  (SVPSF) 

describing  the  effects  of  the  optical  imaging  system. 
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Ideally,  it  is  desirable  that  the  point-spread  function  for 
a Dirac  delta  function  6(x,y),  in  which  case  the  image 
g ( x ,y ) equals  the  object  f(x,y)  in  the  absence  of  noise. 
Furthermore,  if  h(e,n;x,y)  is  a function  only  of  the 
differences  between  respective  coordinates,  that  is 

h(e,n;x,y)  = h(x-e,y-n)  (2.2) 

the  PSF  is  said  to  be  spatially  invariant  or  isoplanar.  In 
all  other  situations,  the  FSF  is  said  to  be  spatially 
variant  or  anisoplanar.  The  physical  meaning  of  a 
spatially  invariant  FSF  is  that  the  blurring  degradation  is 
unchanged  across  the  image  plane  and  the  image  and  object 
are  mathematically  related  via  a two-dimensional 
convolut ion 


g (x,y) 

or  equivalently 

g (x,y) 


h (x-e ,y-n) f (e , n) dcdn+n (x ,y) 


h(e  ,n)  f (x-e ,y-n)dedn+n  (x,y) 


(2.3) 


(2.4) 


These  convolution  integrals  can  be  Fourier  transformed  to 
yield 


G (u , v)  = H (u, v) F (u,v) +N (u, v)  (2.5) 


where  the  capital  letters  denote  the  Fourier  transform  of 


the  respective  function  represented  by  lower  case  letters 
and  (u,v)  is  spatial  frequency.  In  the  discrete  domain 
Eq.  (2.5)  can  be  expressed  by  a discrete  Fourier  transform 
(DFT ) equivalent. 

To  process  image  signals  on  2 digital  computer  we  need 
a spatially  aiscrete  form  of  signal.  Equation  (2.4)  can  be 
represented  as  a aiscrete-discrete  system  [2-1, 2-2)  by  a 
matrix.  This  matrix  can  be  represented  as  a vector  by 
lexicographically  ordering  the  column  of  the  matrix,  i.e., 
the  (i,j)th  element  of  the  M *N  matrix  is  the  l(j-l)m+l]th 
element  of  the  vector.  This  ordering  permits  the  use  of  a 
very  simple  matrix  model 

g = Hf  + n (2.6) 

2 

where  £ is  an  N *1  recorded  or  measured  image  data  vector 
2 

jE  is  an  N *1  original  object  vector 
2 

n is  an  N *1  additive  Gaussian  noise  vector 
2 2 

H is  an  N *N  blurring  matrix  which  is  a 

transformation  matrix  representing  the 
blurring  degradation. 

Thus  the  linear  restoration  matrix  model  is  as  shown  in 
Fig.  2.2.  lhe  additive  nature  of  the  noise  in  Eq.  (2.6)  is 
a model  for  thermal  noise  and  amplifier  electronics  noise 
in  image  sensors.  This  additive  noise  is  often  itself 
modeled  as  a Gaussian  process.  Since  this  model  represents 

the  physical  reality  well  and  is  mathematically  tractable, 
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it  is  the  most  conventional  practical  model. 


2.2.2  Mu  1 1 ipl icat ive  Image  Noise  Model 

A block  aiagram  of  a multiplicative  noise  model 
12-b, 2-7 ,2-fa)  ana  the  associated  restoration  filter  is 
illustrated  in  Fig.  2.3.  Its  mathematical  expression  is 

g(x,y)  = [f (x,y)«h(x,y) ] *n (x,y)  (2.7) 


f 


. 


I 

I 


or 


g (x,y) 


h (c ,n;x,y) 


f (c , n ) drdn 


n (x,y) 


(2.8) 


where  h(x,y)  is  the  FSF  of  the  linear  system,  f(x,y)  and 
g(x,y)  are  the  obejct  and  degraded  image  functions 
respectively.  Here  n(x,y)  denotes  signal- independent 
Gaussian  noise,  and  ® denotes  two-dimensional  convolution. 
T.  Yatagai  (2-9)  has  used  this  model  for  speckle  noise  in 
the  sense  that  standard  deviation  of  the  speckle  is  equal 
to  its  mean. 


2.2.3  Additive  Signal-Modulated  Image  Noise  Model 

The  adaitive  signal-modulated  image  noise  model  and 
the  associated  restoration  filter  is  illustrated  in 
Fig.  2.4.  Its  mathematical  expression  is 

g(x,y)  = f (x,y)»h(x,y)+clf (x,y)»h(x,y) ] *n(x,y) . (2.9) 
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Restoration 

filter 


Block  diagram  of  the  linear  additive 
Gaussian  noise  model  and  the  associated 
restoration  filter 


Blurring 

matrix 


Restoration 
f i Iter 


Block  diagram  of  the  linear  restoration 
vector  model  with  additive  noise 


Restoration 

filter 


Block  diagram  of  the  multiplicative 
model  and  the  associated  restoration 
filter 


Here  C is  generally  a memoryless  point  nonlinear  function, 
n(x,y)  is  the  signal- independent  noise  and  » denotes 
two-dimensional  convolution.  Film-grain  noise  and  magnetic 
tape  recording  noise  are  accurately  modeled  as  additive 
signal-modulated  noise  [2-10,2-11], 


2.2.4  Poisson  Image  Noise  Model 


Because  of  the  quantum  nature  of  light,  photons  arrive 
at  random  times  and  give  rise  to  a fundamental  graininess 
in  detected  images  at  low  light  levels.  The  graininess 
tends  to  obscure  the  detection  of  fine  details  and  faint 
contrasts,  thus,  a large  number  of  photons  is  required  for 
high  quality  imaging.  The  emission  of  photons  is  governed 
by  a Poisson  random  process  [2-19],  hence,  we  label  the 
noise  as  Poisson  noise.  Because  a detailed  mathematical 
model  for  Poisson  noise  is  given  in  section  2.4,  only 
qualitative  comments  on  Poisson  noise  are  given  here  for 
the  completeness  of  this  section. 


Poisson  noise  is  another  basic  type  of 
signal-dependent  image  noise  which  is  quite  different  from 
film-grain  noise  and  speckle  noise.  The  signal-dependent 
nature  of  Poisson  noise  is  associated  with  the  fact  that 
the  variance  of  the  Poisson  probability  distribution  is 
equal  to  its  mean.  If  the  signal  information  received  by 
an  array  of  photo  detector  elements  is  contained  in  the 
mean  number  of  events  recorded  by  each  element,  then  the 


Poisson  distribution  of  these  events  implies  that  Poisson 
noise  is  a form  of  signal-dependent  noise.  All  low-level 
photon  resolved  signals  are  examples  of  signals  corrupted 
by  Poisson  noise.  These  situations  occur  in  scintillation 
camera  imaging,  medical  imaging,  astronomical  imaging,  and 
low  light  level  television  systems. 

2 . 3 Image  Models 

Three  of  the  more  detailed  image  models  are  discussed 
in  this  section.  Better  estimates  of  statistical  images 
should  come  from  more  accurate  image  models. 

2.3.1  Gaussian  Image  Model 

On  the  basis  of  physical  arguments  and  mathematical 
tr actabil ity , Hunt  [2-12,2-13)  developed  one  of  the  most 
accurate  image  models.  The  image  is  modeled  as  a 
multivariate  Gaussian  process  with  nonstationary  ensemble 
mean  and  with  stationary  covariance.  The  image  vector  a 
pr ior i probability  density  function  (pdf)  is 

p(f)  = ((2ir)N/2|Rf  |,,)"1exp{-j(f-I)TRj1  (f-?)  } (2.10) 

where  F=E[f]  is  the  nonstationary  mean  vector 
— — T 

Rf=El  (f-_f) (£“£)  ) is  the  stationary  covariance  matrix 
| R f ) is  the  determinant  of  Rf  and  T denotes 
transpose. 

Lower  case  p denotes  probability  density  function  (pdf)  and 
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capital  case  p denotes  probability  throughout  this  thesis. 
Equation  (2.1U)  describes  a random  process  of  stationary 


I 


» 


fluctuations  about  a nonstationary  mean  value.  Hunt 
derived  this  image  model  based  on  a heuristic  analysis  from 
following  "thought"  experiment.  Suppose  that  several 
thousands  of  photographs  with  similar  statistical 
properties  (such  as  driver's  license  photographs)  were  used 
to  calculate  an  ensemble  mean  image.  Each  face  is 
positioned  in  approximately  the  same  way  in  each  image 
frame.  Clearly  such  an  ensemble  mean  image  would  not 
consist  of  a uniform  shade  of  gray  indicating  spatial 
stationary.  More  likely,  the  mean  image  probably  consists 
of  an  oval  region  where  the  face  is  expected  to  be  and  some 
dark  spot  where  the  eyes,  nose,  and  other  facial  features 
are  expected  to  be.  Thus,  images  are  generally 
nonstationary  in  first  order  statistics  and  are  described 
by  a spatially  non-stat ionary  mean  vector  £.  The  ensemble 
mean  is  strongly  dependent  upon  the  context  which  is 
established  by  the  sample  mean  of  the  image  to  be  modeled. 
The  ensemble  is  called  a context-dependent  image  ensemble. 
Wintz  12-16]  has  shown  that  images  may  have  identical 
covariance  statistics  and  the  same  constant  mean  intensity, 
but  be  completely  unrelated  in  context.  Therefore , the 
context-dependent  ensemble  properties  are  portrayed  most 
strongly  in  the  mean  vector  £,  since  this  vector  has  the 
gross  structure  that  represents  the  context  of  the  ensemble 
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from  which  the  sample  vector  jE  is  drawn. 


I 


The  spatial  statistics  f and  of  Eq.  (2.12)  do 
not  represent  a "second-oraer " stationary  random  process  if 
the  context-dependent  ensemble  mean  is  assumed,  since  f is 
assumed  to  have  nonstationary  gross  structure  that  depends 
upon  the  ensemble  and  its  context.  The  covariance 
statistics  can  be  described  as  spatially  stationary 
fluctuations  about  a spatially  nonstationary  mean  vector 
jf.  The  random  process  associated  with  the  image  ensemble 
is  not  ergodic  in  the  mean,  since  the  ensemble  average  of 
the  context-dependent  ensemble  is  not  equal  to  the  spatial 


average  of  an 

ensemble  member. 

However,  the  process 

can 

have  a stationary  autocovariance. 

If  an  image  f(j,k)  can  be 

described  as  being 

the 

sum 

of  two  components,  a low  frequency  or  blurred 

component 

f ( j , k ) and 

a high-frequency 

component  s(j,k) 

of 

the 

fluctuations 

about  f,  i.e. 

f ( j , k ) = f(j,k)  + s(j,k) 

(2. 

ID 

then  F(  j ,k) 

represents  the 

nonstationary  mean 

and 

the 

variance  of 

the  difference 

image  between  the  image  f 

and 

the  nonstationary  mean  image  f 

is  approximately 

equal 

to 

that  of  the 

high-frequency 

component  s( j ,k)  . 

Thus , 

the 

ensemble  mean 

of  the  random 

process  for  an 

image 

is 

nonstationary  and  carries  the  low  spatial  frequency  gross 
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features  of  an  image,  while  the  covariance  properties  about 


this  ensemble  mean  represent  random  perturbations  carrying 
the  detailea  image  structure.  This  image  model  is  one  of 
the  most  sophisticated  and  will  be  used  in  deriving  optimum 
nonlinear  filters  in  later  chapters. 


.acian  Image  Model 


Trussel  et  al.  [2-13]  using  a nonlinear  least-square 
fit  technique  found  that  a Laplacian  pdf  had  a better  fit 
to  experimental  image  pdf's  than  the  Gaussian  pdf.  Thus 
the  model  of  Eq.  (2.9)  can  be  restated  as 


PL(f)  = k2exp{-/2  [ (f-f)TRf  (f-f)P) 


(2.12) 


where 


f = E[f] 

Rf  = E[  (f-f)  ( f-f ) T] . 

This  process  is  similar  to  the  Gaussian  model  because  it 
contains  a nonstationary  mean  and  stationary  covariance. 
The  square  root  of  the  term  in  brackets  in  Eq.  (2.12)  makes 
the  model  Laplacian  instead  of  Gaussian.  The  MAP  filter 
and  other  restoration  methods  are  later  derived  for  both 
the  Gaussian  pdf  and  Laplacian  pdf. 


2.3.3  Lebeoev 1 s Composite  Image  Model 


This  image  model  repesents  a completely  different 
approach  from  Hunt's  image  model.  Lebedev  and  Mirkin 
called  their  moael  a "composite  model  of  an  imago  fragment" 


which  includes  the  nonstationary  statistical  properties  of 
an  image  [2-14,2-151.  They  model  an  ensemble  of  images  as 
a random  field  with  an  n-dimensional  point  pdf  p(f)  where 
- = l f i' f 2 ' ' • * • f * They  decompose  the  image  statistically 
into  M classes  of  fragments,  whose  structure  is 
distinguished  by  the  type  of  correlational  links  between 
pixels.  Some  classes  are  formed  by  fragments  with 
isotropic  structure;  others  are  found  by  fragments  with 
some  anisotropy. 

Let  p(f)  be  the  pdf  of  a fragment  of  image  f,  on  the 
condition  that  the  fragment  belongs  to  the  class  9 
(6  = 1,2,...  ,M)  . Denoting  the  a pr ior i density  of  the 
classes  by  it(0),  wc  have 


M 


and 

M 

p(f)  = £ rr(0)p0(f).  (2.14) 

0 = 1 

Expression  (2.14)  is  a matrix  density  decomposition  of  p(f) 
in  terms  of  pQ(f) , 0 = 1,2,... ,M  [2-18].  This 

representation  is  especially  useful  when  p(f)  is 
approximated  closely  by  a small  collection  of  Gaussian 
d istr ibut ions 
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(2.15) 


p0(f)  = | R0  r1exp{-ifTRQ1f  } 


where  Rfl  is  the  covariance  matrix  corresponding  to  the 
class  0 fragment  image. 

Using  this  image  model  for  developing  a spatial 
restoration  filter  with  a maximum  likelihood  (ML),  maximum 
a poster  ior  i (MAP)  or  Wiener  criterion  leads  to  a 
multicategory  filtering  problem  because  of  the 
decomposition  properties  of  the  image  model.  This 
composite  image  technique  seems  to  be  a good  model  for  thc- 
local  nonhomogeneous  information  in  the  image  signals, 
hence  an  optimal  filter  can  be  a local  adaptive  filter. 
Although  this  model  is  not  used  in  the  results  presented, 
it  is  believed  that  using  this  model  with  the  MAP  criterion 
may  yield  good  results  in  future  work. 

2.4  Inhomogeneous  Poisson  Process  Photon  Counting  System 
Model 

2.4.1  Photon  Counting  System 

In  many  practical  situations  a detected  image  can  be 
modeled  as  a photon  counting  system  illustrated  in  Fig.  2.5 
with  its  corresponding  block  diagram  shown  in  Fig.  2.6. 
For  ease  of  notation,  we  use  a lexicographic  ordering 
vector  notation  in  which 


H 


(Aberrations, 
blur  diffraction 
turbulence) 


Estimate  f 


Figure  2.5.  Two-dimensional  photon  counting  system 


we  first  derive  the  statistics  of  a single  counter  and 
later  extend  it  to  a large  array. 

2.4.2  Single  Photon  Counter 

From  Fig.  2.5  we  first  assume  the  H matrix  is  the 
identity  matrix,  and  assume  only  one  counter  instead  of  a 
vector  array  as  in  Fig.  2.7.  According  to  the 
semi-classical  theory  of  photon  detection  [2-19,2-26],  the 
probability  that  g^  photon  events  occur  for  a given  fixed 
intensity  f^  is 

g -Af. 

(Af.)  e 

p ] f ) = i , (2.16) 

g i 

i 

where  A is  a constant  rate  parameter.  By  direct  summation, 
the  conditional  mean  ana  variance  of  g^^  for  constant  f ^ is 

qA  = Xfi#  (2.17) 


and 


(2.18) 


- X = Average  # of  photon  counts 
f^  Intensity  unit 

With  the  low  light  level  image  signals  in  which 
interested,  we  have  12-21,2-26] 


WT  > > 1 , 


or  equivalently, 

g^  < < WT , (2.21) 

where  W is  the  temporal  bandwidth  and  T is  the  integration 
time,  i.e.  the  mean  number  of  counted  photo  events  is 
small  compared  to  the  time-bandwidth  product  of  the  light. 
This  condition  is  always  satisfied  for  natural  thermal 
sources  encountered  in  practice.  In  such  cases,  Mandel 
[2-26]  has  shown  that  the  count  fluctuations  are 
predominantly  Poisson  shot  noise  due  to  the  discrete  nature 
of  wave-detector  interaction,  rather  than  classical 
"fluctuation  noise"  associated  with  the  random  nature  of 
the  image  intensity.  Thus  the  count  registered  on  the  ith 
counter  is  a Poisson  random  variable  with  mean  g ^ and  pdf 
expressed  by  Eq.  (2.16). 

From  the  linear  transformation  d^=ag^  with  a a 
constant  display  scale  factor,  we  have 

^i  = »gi  = a A f ^ , (2.22) 


and  p(d.|f.)  is  given  by 


I 

1 


the  processed  display  image  signal  the  same  as  the  observed 
noisy  image. 

Before  we  aefine  an  r.m.s.  signal- to-noise  ratio,  it 

is  useful  to  discuss  some  details  about  signal-to-noise 

ratios  (SNR's)  in  general.  A study  of  SNR's  at  different 

points  in  a system  enables  us  to  pinpoint  the  significant 

contributions  to  the  noise.  It  is  also  a simple  criterion 

for  the  design  of  systems  to  minimize  from  the  noise 

degradation  and  thus  it  provides  a measure  of  the 

"noisiness"  of  a system.  In  most  cases,  the  SNR  criterion 

is  applied  with  signal- independent  zero  mean  additive 

noise.  If  the  signal  and  noise  are  dependent,  then  the  SNR 

is  difficult  to  define  because  the  cross  correlation  and 

other  moments  are  nonzero.  Poisson  noise  is  a case  of 

signal-dependent  noise  because  the  variance  in  Eqs.  (2.18) 

and  (2.23)  depends  on  the  signal.  In  order  to  compare  the 

noisiness  of  images  with  Poisson  noise,  we  will  define  an 

r.m.s.  SNR  denoted  by  (SNR)rms.  Because  an  image  signal  is 

a random  process  in  space  and  time,  we  must  define  an 

ensemble  mean  (SNR)  as  the  ensemble  average  of  the 

rms 

(SNR)  over  the  random  image  field. 

rms 

In  the  case  of  Poisson  noise,  the  (SNR)rms  is 


(2.25) 
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and  the  ensemble  mean  (SNR) 


is 


rms 


ISNR,rms  * ^ * X£i 


(2.26) 


Thus,  (SNft) is  proportional  to  the  square  root  of  the 
signal,  and  is  signal-dependent. 


This  behavior  is  quite  different  from  additive  noise, 
multiplicative  noise,  or  film-grain  noise.  Figure  2.8  is  a 
digital  simulation  of  a low-level  one-dimensional  image 
signal  whose  peak  SNR  is  approximately  18.  It  is  clear 
from  the  illustration  of  Fig.  2.8,  that  at  low  signal 
levels  the  noise  is  statistically  nonstationary  and 
non-Gauss ian.  As  the  signal  becomes  photon  resolved  at 
very  low  intensities,  there  is  little  resemblance  to  the 
classical  signal.  However,  at  higher  signal  levels,  the 
noise  becomes  more  Gaussian  [2-27]. 

2.4.3  Statistics  of  Array  Counters 


An  array  counter  model  for  non-blurred  image  signals 
with  Poisson  noise  is  shown  in  Fig.  2.9.  For  one  single 
counter  the  conditional  density  is  given  by  Eq.  (2.16). 
For  an  array  of  counters,  we  must  find  the  joint  ensemble 
statistics  for  a given  object  vector  f^.  Some  assumptions 
are  necessary  to  find  these  joint  ensemble  Poisson 
statistics.  Walkup  [2-20,2-21],  Clark  [2-22],  and  Wang 
[2-23]  have  shown  that  given  f,  the  joint  ensemble  photon 


Parameter  X 


Display 

intensity 
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counting  statistics  at  the  N detectors  of  Fig.  2.5  are 
independent  Poisson  variates  whose  mean  is  proportional  to 
the  mean  intensities  at  each  of  the  N disjoint  detectors. 
For  this  to  be  true,  the  time-bandwidth  product  WT  must 
satisfy  WT>>1  where  W is  the  temporal  bandwidth  of  the 
optical  image  signal  and  T is  integration  time.  This  is 
the  usual  case  in  this  thesis.  From  the  WT  >>1  assumption, 
Mandel  [ 2— 2 fo ] also  has  shown  that  the  counts  registered  by 
the  N detector/counters  may  be  taken  to  be  statistically 
independent,  since  classical  fluctuation  noise  is 
negligible  for  light  with  a low  degeneracy  parameter,  when 
Poisson  shot  noise  predominates  in  the  photon  count 
fluctuations.  The  degeneracy  parameter  describes  the 


variance  characteristics  of  the  photon  counts,  and  it  is 
defined  as  the  ratio  of  classical  fluctuation  noise  to 
Poisson  shot  noise.  Thus,  all  g^  are  independent  for  a 
given  t (i.e.  every  Poisson  generator  is  independent)  and 
each  g^  depends  only  on  its  corresponding  f ^ . Also,  we 
assume  that  individual  detectors  have  a smaller  scale  than 
the  spatial  intensity  variations  of  the  image  so  that  no 
loss  of  information  results  from  the  sampling.  We  also 
initially  assume  that  background  intensities  and  thermal 
noise  in  the  measurement  system  is  negligible  compared  to 
the  Poisson  noise.  Representing  the  array  values  by 

2=  [g1,g2, . . . ,gNJT  and  f = [f^  f2  , . . . , fN)T  (2.27) 


i 


we  have 


p(£|f)  = p(g1lf)p(g2lf) . . .p(gNlf)  (2.28) 


Now,  each  gi  depends  only  on  its  corresponding  f^,  thus 


g.  -Xf . 
(Xf.)  xe  x 

p(2|f)  = n — 


g . ! 
1 


(2.29) 


Now  from  Eqs . ( 2 . 22 ) - (2 . 24 ) we  have 


p(d|f)  = ptdj^  |f)p(d2|f)  . . .p(dN|f) 

= p(d1|f1)p(d2|f2)...p(dN|fN) 


(2.30) 


p (d  | f ) = n 
i 


-Xf. 

(Xfi)  a e 1 

T. 

a( — -)  ! 
a 


(2.31) 


From  this  conditional  density  an  MAP  estimate  is  derived  in 
Chapter  4. 


2.4.4  The  Quantum  Limitations  of  Photon  Resolved  image 

F.i  ~nals 

The  information  content  of  a finite  amount  of  light  is 
limited  by  the  finite  number  of  photons,  by  the  random 
character  of  their  distribution,  and  by  the  need  to  avoid 
false  alarms.  These  limitations  mean  that  a considerable 
number  of  photons  is  required  to  delineate  the  fine  detail 
of  images.  Low  light  level  image  signals  conspicuously 
suffer  from  these  quantum  limitations  (2-24). 
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(i)  Discreteness  of  Light  Quanta 


The  number  of  picture  elements  N required  for  a well 

resolved  image  signal  often  lies  in  the  range  of  10^-107 

8 9 

[2-24].  Hence,  10  -lu  photons  are  needed  to  delineate  the 
location  and  brightness  of  pixels,  assuming  that  100 
photons  per  pixel  are  arranged  in  a precise  array. 
Unfortunately,  nature  does  not  work  in  so  orderly  a 
fashion,  and  photons  arrive  at  random  times  and  places  and 
give  rise  to  a fundamental  graininess  in  any  detected 
image.  This  grain  noise  obscures  the  detection  of  fine 
detail  in  low  contrast  images  and  is  called  Poisson  noise 
at  these  low  light  levels. 

(ii)  Random  Character  of  the  Photon  Distribution 

Suppose  that  the  average  number  of  photons  arriving  in 

a given  area  is  ng , and  that  a number  n^ ,n2 » . . . »n^  are 

distributed  around  ng  in  such  a fashion  that  the  average 

value  of  (n^-ng)  is  also  ng.  Then  the  mean  of  the  photon 

1/2 

count  is  the  same  as  the  variance  and  the  (SNR)rmsis  ng  . 

If  we  need  to  detect  a 1%  contrast  variation  in  a 
signal,  we  require  that  the  noise  level  defined  as  the 
r.m.s.  deviation  of  the  mean  number  of  photons  also  not 

exceed  1%.  This  can  be  achieved  by  having  an  average  of 

4 1 

10  photons  falling  on  each  pixel  of  the  image.  The  r.m.s.  I 

4 2 ! 

deviation  (noise)  is  then  the  square  root  of  10  , or  10  , I 


r.-’x^TJ-v 
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ana  the  ratio  of  this  random  deviation  to  the  average  will 
-2 

be  1U  , or  1%.  In  general,  if  we  want  to  detect  an  image 

signal  with  contrast  C we  must  increase  the  number  of 

photon  counts  proportional  to  In  addition,  we  must 

C 

guard  against  false  alarms,  that  is,  the  mistaking  of  any 

■ 

particularly  random  fluctuation  for  the  real  signal  to  be 
detected . 

(iii)  False  Alarms 


False  alarms  and  spurious  visual  patterns  may  arise 
from  the  random  character  of  the  photon  distribution  and 
not  from  the  original  scene  itself.  If  we  define  the  mean 
photon  count  as  the  signal,  and  the  r.m.s.  deviation  as 
noise,  then  we  can  use  detection  theory  to  find  the  false 
alarm  probability.  Figure  2.10  shows  the  distribution  of 
noise  fluctuations  around  the  mean  value  of  a parameter. 
The  ordinate  is  the  probability  density  and  the  abscissa  k 
is  plotted  in  un'ts  of  the  r.m.s.  deviation.  The  second 
abscissa  scale  n is  a particular  example  for  which  the 
average  number  of  photons  is  900  and  its  r.m.s.  deviation 
is  30.  From  Fig.  2.10,  we  can  calculate  the  false  alarm 
probability  given  a signal  which  is  k units  above  its  mean 
value.  Table  2.1  gives  the  probability  that  noise 
fluctuations  will  exceed  the  mean  value  of  the  background 
by  an  integer  number  of  units  of  the  r.m.s.  value  of  the 
noise . 
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Figure  2.10.  Probability  distribution  of  a RMS 

deviation  of  photon  counts  about  its  mean 


Table  2 . 1 

False  alarm  probability  of  exceeding  various 
values  of  k 


k 


prob.  of  exceeding  k 


If  we  locate  the  signal  at  k=5,  then  we  find  that  only 

15%  of  the  time  will  the  signal  appear  k<4 . It  will  exceed 

k=4  85%  of  the  time  on  the  average  and  be  judged  a real 

signal.  Hence,  we  usually  choose  k=5  to  avoid  false  alarms 

in  order  to  give  a reasonable  reliability  to  our 

observation.  The  ratio  of  the  r.m.s.  deviation  to  the 

average  background  brightness  varies  as  n0  /n0=l/n0  , where 

nQ  is  the  average  number  of  photons  in  the  background. 

Hence,  it  will  be  necessary  to  increase  n by  a factor  of 

k in  order  to  decrease  the  ratio  l/nQ  by  k.  In  summary, 

2 

the  density  of  photons  required  varies  as  k in  order  to 
avoid  false  alarms.  In  general,  the  expression  for  the 
total  number  of  photons  required  to  detect  a contrast  C 

A 

where  C=  B/B  and  0<C£l  is 

1 2 

Total  number  of  photons  = N — =■  k (2.32) 

C 

Here,  N is  the  total  number  of  pixels  in  the  picture  and 

reflects  the  discrete  nature  of  the  photons.  The  factor  -K 

C 

is  a consequence  of  contrast  the  C and  the  random  character 

2 

of  photon  distributions,  the  factor  k reflects  both  the 
random  character  of  photon  distribution  and  the  need  to 
avoid  false  alarms.  The  expression  (2.32)  is  only  for  the 
case  in  which  we  do  not  have  any  a priori  information  about 
the  image.  However,  we  usually  assume  that  the  image  is  a 
Markovian  random  field  with  correlation  coefficient  p in 
image  restoration  work.  Thus,  we  should  be  able  to  obtain 
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better  results  than  obtained  in  Eq.  (2.32) 


■ 

I 


2.4.5  Comparison  Between  Linear  Add i t i ve  Gaussian  Noise 
Degraded  Image  Signals  and  Poisson  Noise  Degraded  Image 
Signal s 

Figure  2.11  illustrates  experimental  results  of 
pictures  with  Poisson  noise  and  Gaussian  noise  for 
comparable  (SNR) 

rms 

The  upper  left  image  (A)  in  Fig.  2.11  is  the  original 

ideal  object.  The  upper  right  image  (B)  contains  Poisson 

noise  with  some  constant  amplification  and  (SNR)  „=6  db. 

rms 

The  lower  left  image  (C)  has  additive  Gaussian  noise  with 

(SNR)  also  6 db.  The  lower  right  one  (D)  also  has 
i ms 

linear  additive  Gaussian  noise  with  (SNR)  approximately 

rms 

1 ob.  From  this  experimental  result,  we  have  demonstrated 
that  images  with  Poisson  noise  are  more  severely  degraded 
than  images  with  linear  additive  Gaussian  noise  for 
comparable  (SNR)  even  though  the  (SNRl  „ is  b db  lower 
than  that  of  an  image  with  Poisson  noise.  The  Poisson 
noise  has  obliterated  almost  completely  the  detailed  edge 
and  contrast  of  the  face.  The  image  takes  on  a mottled 
appearance  which  depends  on  brightness  whereas  the  image 
degraded  by  Gaussian  noise  appears  uniformly  degraded  with 
some  edge  contrast  still  discernible  at  comparable  (SNR). 


4 


l 
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2 . b 


Conclusion  s 


We  have  presented  some  system  models  for  image  noise 
and  three  image  models.  We  have  discussed  the  difference 
between  the  Poisson  image  noise  model  and  other  important 
image  noise  models  from  a system  point  of  view.  The  model 
which  we  have  developed  in  the  section  2.4  is  used  in  later 
chapters  for  developing  and  implementing  an  MAP  filter  for 
the  Poisson  noise  model.  Appendix  A describes  a computer 
algorithm  which  is  used  to  generate  Poisson  noise  for  all 
simulations  throughout  this  thesis.  The  quantum 
limitations  of  photon  resolved  image  signals  are  presented 
in  order  to  further  understand  the  physical  meaning  and 
causes  of  Poisson  noise  in  low  light  level  image  signals. 
By  modeling  the  low  level  signals  as  an  inhomogeneous 
Poisson  process,  a very  accurate,  complete  model  for  many 
physical  systems  including  medical  imaging,  astronomical 
imaging  is  developed. 
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CHAPTER  3 


REVIEW  OF  LINEAR  AND  NONLINEAR  OPTIMAL  IMAGE  RESTORATION 
FILTERS  FOR  THE  IMAGE  NOISE  MODELS 

3 . 1 Introduction 

In  this  chapter,  we  review  previous  work  on  linear  and 
nonlinear  restoration  filters  with  Poisson  image  noise 
models  and  other  important  image  noise  models.  These  image 
spatial  restoration  fillers  are  developed  based  upon  the 
system  models  for  image  noise  given  in  Chaper  2.  The 
relationship  of  these  fillers  to  the  work  presented  in 
later  chaptes  is  also  discussed. 

3 . 2 Spatial  Restoration  Filters  for  the  S ignal-I ndependent 
Noise  Model  and  S igna 1 -Dependent  Noise  Model 

3.2.1  Spatial  Restoration  Filter  for  Signal-Independent 
Noise Linear  Additive  Gaussian  Noise  Model 

Most  previous  work  in  image  restoration  is  based  on 
the  model  of  section  2.2.1.  Different  filters  such  as  the 
inverse  filter,  constrained  least-squares  filter, 
parametric  Wiener  filter,  homomorphic  filter,  maximum 
entropy  filter,  and  pseudo- inver se  filter  have  been 
developed  under  various  criteria  [3-1, 3-2].  Hunt  and 


Trussell  [3-3,3-4,3-5]  have  also  developed  a nonlinear  MAP 
filter  based  on  the  direct  maximization  of  the  posterior 
density  function  with  a linear  additive  image  noise  model. 


3.2.2  Wiener  Spatial  Filters  for  Signal-Dependent  Noise 

Model Mu 1 tipi icative  Noise  , Additive  Signal -Modulated 

Noise 


The  models  in  sections  2.2.2  and  2.2.3  are 

signal-depdnent  noise  models.  Walkup  et  al . [3-7],  Hondo 

et  al . [ 3—8  J and  Yatagai  [3-9]  developed  an  optimal  spatial 

filter  in  the  sense  of  minimizing  the  mean  of  the  squared 
~ 2 

error  [ f (x  ,y ) -f ( x ,y) ] in  the  manner  of  a Wiener  filter. 
From  the  minimum  mean-square  error  criterion,  the 
orthogonality  principle  is  developed,  leading  to  a spatial 
Wiener  filter  w=Rfg(Rggl  where  Rfg  and  Rga  are  the 

cross-covariance  matrices  between  object  and  image  and  the 
cova'  xance  matrices  of  the  image  respectively  [3-12].  If 
g(x,y)  and  f(x,y)  are  spatially  wide-sense  stationary 
random  processes,  then  the  two-dimensional  Wiener  filter 
has  a spatial  frequency  domain  transfer  function 


uj(u,v)  = 


♦fq(U'V) 


(3.1) 


where  *fg  and  *gg  represent  the  cross-spectral  densities  of 
the  image  g(x,y)  with  the  object  f(x,y)  and  the  spectral 
density  of  the  image  g(x,y),  respectively. 


Naderi  and  Sawchuk  [3-10]  also  developed  a better 
spatial  adaptive  Wiener  filter  for  signal-dependent 
film-grain  noise  using  a more  accurate,  complex  noise 
model.  Walkup  £t  al . [3-11]  developed  an  MAP  spatial 
filter  for  signal-dependent  noise  models  such  as  film-grain 
noise  and  magnetic  recording  noise.  Their  filter  involved 
scalar  rather  than  vector  processing  and  did  not  include 
blurring  degradations. 

3 . 3 Linear  , Invariant  Least-Squares  Restoration  Filter  for 
the  Poisson  Image  Noise  Model 

Goodman  and  Belsher  [3-13]  have  modeled  a low  light 
level  imaging  system  with  blurring  and  Poisson  noise  by  the 
continuous  system  shown  in  Fig.  3.1.  This  system  is  a 
special  case  of  the  discrete  restoration  model  discussed  in 
section  2.4.  The  detected  data  is  represented  as 

N 

g(x,y>  - E «<x-x  ,y-y  > <3.2> 

n=l 

where  6 ( - , * ) is  a two-dimensional  Dirac  delta  function 
(xn»yn)  represents  the  location  of  the  nth  photo 
event 

N is  the  total  number  of  photo  events  produced  by 
the  blurred  image  b(x,y) 

In  the  expression  (3.2),  N,  xn  and  yn  are  all  regarded  as 
random  variables. 
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From  the  semi-classical  theory  of  photo  detection 


[3-14,3-21],  the  probability  that  N events  occur  in  an  area 
A on  the  detector  is  taken  to  be  Poisson, 


Pa(N)  = 


_ ’ 

, , X (x, 
L A 


y ) dxdy 


>1-  II  A (x, 

1 A 


y ) dxdy 


(3.3) 


where  A(x,y)  is  a rate  function 


i(x.y>  - nMiLil  T 
h v 


n is  the  quantum  efficiency 
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h is  plank's  constant  = 6.624*10  W-sec/HZ 
v is  the  mean  optical  frequency 
t is  the  detector  integration  time 


:,y ) = j,  h(c,n;x,y)f  (e, 


n)dedn 


(3.4) 


(3.5) 


Because  the  image  f(x,y)  is  a random  process,  the  rate 
function  A(x,y)  is  a random  process.  Thus  the  whole 
process  is  called  a compound  Poisson  process  or  doubly 
Poisson  process  [3-14,3-15].  The  event  locations  (xn»yn) 
are  independent  variables  for  different  n's  for  Mx,y) 
given,  and  it  has  pdf. 


x(vV 


II  A(X, 


(3.6) 


y)dxdy 


In  Fig.  3.2,  we  illustrate  a one-dimensional  typical  object 
intensity  distribution  and  a corresponding  typical  detected 
image. 


Based  on  the  model  of  Eq.  (3.2)  and  Fig.  3.1,  Goodman 
and  Belsher  [3-13]  developed  an  linear,  space- invar iant , 
least-square  restoration  filter  ( LMMSE  filter).  The  filter 
is  derived  in  the  Fourier  transform  domain  using  the 
orthogonality  principle.  The  form  of  the  filter  is 

N V (u,v)$f(u,v) 

W(u,v)  = — — (3.7) 

1+N  I ^ ( u , v)  I (u,v) 

where  % is  the  two-dimensional  Fourier  transform  of  h(x,y) 
4>f  is  the  spectral  density  of  the  object  f(x,y) 

N is  the  mean  number  of  photon  counts 

* denotes  complex  conjugate 

(u,v)  are  spatial  frequency  variables 

Equation  (3.7)  can  be  rewritten  as 


W (u,v) 


v (u,v) 
lv (u,v) | 2+i 


(3.8) 


where  a=N4>^(u,v).  This  filter  can  be  implemented  in  the 
Fourier  domain  using  the  FFT.  Goodman  and  Belsher  did  not 
apply  this  filter  to  two-dimensional  picture  data,  so  it  is 
implemented  in  Chapter  7 and  compared  with  the  MAP  filter 
described  in  Chapter  4 and  b.  From  Eq.  (3.8),  if  a is  very 


large,  then  W(u,v)=H  *( u,v ) which  is  the  ideal  inverse 
filter  in  the  absence  of  noise.  Hence,  the  LMMSE  filter 
should  work  very  well  at  high  values  of  a.  However,  for 
the  low  light  levels  in  which  Poisson  noise  dominates,  a is 
very  small.  Therefore,  the  performance  of  the  LMMSE  filter 
is  expected  to  be  poor  due  to  ill-conditioning  in  the 
deconvolution  process.  Furthermore,  the  LMMSE  filter 
assumes  that  the  object  signal  and  noise  must  be  at  least 
wide-sense  stationary  and  requires  knowledge  of  the  blur 
function  and  object  and  convariance  functions.  It  cannot 
be  extended  in  the  Fourier  domain  to  space-var iant 
blurring.  Because  it  assumes  that  signal  and  noise  are 
stationary,  it  leads  to  a filter  that  tends  to  smooth  edges 
because  of  its  insensitivity  to  abrupt  changes. 

3.4  Nonl  inear  Filter inq  with  the  Square-Root  Transform  for 
the  Poisson  Image  Noise  Model 

Inouye  [3-16]  developed  an  ad  hoc  nonlinear  filter 
with  a square-root  transformation  to  suppress  Poisson  noise 
in  nuclear  medicine  images.  He  assumed  that  the  observed 
data  g(x,y)  is  the  summation  of  the  object  function  f(x,y) 
and  quantum  fluctuation  n(x,y)  as  expressed  by 

g(x,y)  = f (x,y)+n (x,y)  (3.9) 


where  n(x,y)  depends  on  f(x,y)  according  to 


where  a is  a constant  of  proportionality.  Taking  the 


square  root  of  both  sides  of  Eq.  (3.9)  we  get 


IgU/Y)]1*  - (f  (x,y)+n  (x,y)  l*5, 


Igtx.yn11  - . 


(3.12) 


i 


l 


Taking  the  first  two  terms  of  a Taylor  series  expansion  of 
the  above  and  assuming  f ( x ,y ) >>n ( x ,y)  yields  an 
approximation 


(g(x,y)l'i  = f'>(x,y,(l+^iT)  . 

We  now  rewrite  Eq.  (3.13)  as  follows: 

[g(x,y)]i5  = f(x,y)  + n (x,y)  , 


where 


f (x,y)  = f*5  (xf  y)  , 
n(x,y)  = jn (x,y) 
From  Eqs . (3. 10)  and  (3.11),  we  have 


n (x,y ) = 0, 


- - : 


-7r-rr 


(3.13) 


(3.14) 


(3.15) 


I «• 
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and  the  noise  statistics  of  the  square-root  transform  are 
stationary  [3-17].  Thus  this  transform  of  g(x,y)  gives  a 

k 

signal-independent  fluctuating  component  n ( x ,y ) /2f  ( x ,y) 

k 

around  an  average  value  (f(x,y))  . Equation  (3.13)  becomes 

an  approximately  linear  additive  noise  model  and  the  usual 
linear  filtering  technique  is  applied.  Inouye  did  not 
include  any  burring  effects  of  the  imaging  system  and  it  is 
difficult  to  judge  the  level  of  improvement  from  the  line 
printer  pictures  in  his  paper.  In  addition,  this 
square-root  transform  works  only  at  higher  (SNR)  image 
signals,  because  he  assumes  f ( x ,y ) >>n ( x ,y ) . 

3.5  Nonlinear  Optimal  Filter  for  the  Poisson  Image  Noise 
Model 

From  Bayes'  law  we  have  the  posteriori  conditional 
density 

p(d|f)p(f) 

P(f|d)  = , (3.17) 

p (d) 

where  £ is  the  original  image  (object)  to  be  estimated  and 
d is  the  observed  data.  The  use  of  the  posterior  density 
for  estimation  is  well  known  [3-18].  The  minimum 
mean-square  error  estimate  (MMSE)  is  the  mean  of  the 
posterior  density  given  by  E[fjd],  the  maximum  a posterior i 


47 


estimate  (MAP)  is  the  mode  of  the  posteriori  density  as 


expressed  by  max  p(^|d).  The  maximum  likelihood  estimate 
f 

(ML)  given  by  max  p(d|fj  can  be  viewed  as  a special  case  of 
the  MAP  estimate  when  the  a posteriori  density  is  equal  to 
a priori  density  (i.e.  when  p(^f)  has  a uniform 
distribution).  The  MMSE , MAP,  and  ML  estimates  are 
generally  nonlinear,  depending  on  the  form  of  the 
probability  density  functions.  The  MMSE  estimate  also 
needs  the  density  of  the  observed  data  p(d) , but  this  is 
often  difficult  to  obtain  in  practice.  Thus,  a linear 
minimum  mean-square  error  (LMMSE)  estimate  is  commonly  used 
as  described  in  section  3.3.  The  MAP  estimate  tries  to 
find  the  value  of  object  f^  which  maximizes  the  posterior 
density  p(f  |d) . Thus,  it  does  not  need  p(d)  at  all  but 
does  need  p(d  [f)  and  p(f) , while  the  ML  estimate  only  needs 
the  a priori  density  p(d|fj. 

Burke  [3-19]  has  developed  a ML  spatial  filter  based 
on  the  Poisson  noise  model  in  Fig.  3.1.  Its  ML  estimate  is 
obtained  recursively  according  to  the  iteration 


f (n+1 ) _ f(n) 
rk  “ rk 


(3.18) 


where  fk  is  the  kth  component  of  the  object  to  be  estimated 


g . is  the  observed  data  £ of  the  jth  component 
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Hjj  is  ijth  component  of  the  blurring  matrix  H 
3 is  the  quantum  efficiency  constant. 


ans  the  superscript  denotes  the  iterative  step. 

Burke's  experimental  results  are  very  impressive 
compared  to  the  Wiener  filter  results.  Unfortunately,  only 
simple  images  of  very  small  size  (32x32  pixels)  were 
processed  and  the  recursive  Eq.  (3.18)  converges  very 
slowly.  The  ML  estimate  assumes  that  no  information  can  be 
extracted  from  the  a priori  term  p(f) . Modeling  the  object 
as  a nonstationary  random  field  with  a probability  density 
p(f)  to  perform  an  MAP  estimate  can  be  thought  of  an 
extension  of  ML  estimation. 


I I 


i 


4 


3 . 6 Conclusions 

It  is  well  known  that  linear  least  squares  image 
restoration  is  not  optimal  in  the  sense  of  maximum 
likelihood  or  maximum  a posteriori  probability  when  the 
image  statistics  are  Poisson.  Rather,  nonlinear  filtering 
is  required  for  true  optimality.  In  addition,  Fourier 
techniques  cannot  treat  space-variant  imaging.  Thus,  it 
seems  reasonable  that  a nonlinear  filters  should  perform 
better  than  linear  space- invar iant  filters  in  the  presence 
of  signal-dependent  Poisson  noise.  In  subsequent  chapters, 
we  try  to  formulate  and  implement  MAP  filters  in  order  to 
perform  nonlinear  MAP  filter  for  the  low  SNR's  image 
signals  with  Poisson  noise. 
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CHAPTER  4 


RESTORATION  OF  IMAGE  SIGNALS  WITH  POISSON  NOISE 

4 . 1 Introduction 

In  this  chapter,  we  develop  an  optimal  MAP  filter  for 
non-blurred  image  signals  with  Poisson  noise  based  on  the 
model  developed  in  previous  chapters.  There  are  many 
practical  situations  in  which  the  degradation  due  to 
blurring  is  negligible  or  zero.  More  importantly,  we 
develop  a framework  for  MAP  estimation  with  a Poisson  noise 
model  to  evaluate  the  concept  for  future  application  to 
more  complex  problems.  In  section  4.2  we  present  the 
formulation  and  solution  to  the  MAP  estimate  equations.  In 
section  4.3,  4.4,  and  4.5  we  address  the  implementation  of 

an  MAP  estimate  with  different  a pr ior i knowledge.  In 
section  4 . t> , we  describe  a recursive  MAP  filter  for  the 
Poisson  noise  model.  In  section  4.7,  we  discuss  a local 
aoaptive  MAP  filter  and  finally  some  conclusions  are 
presented . 

4 . 2 MAP  Est imate  Formulation  [ 4-10 ] 

Previously  we  derived  a conditional  density  p(d 
displayea  image  data  a given  an  ideal  object  f 
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If)  for 
with  a 


Poisson  noise  mooel.  Also  from  chapter  2,  we  have  the  pdf 
of  the  object  expressea  as  Eq.  (2.1U). 


?rom  the  MAP  estimate  definition,  we  have 

p(d| f)p(f) 

Max  p ( f | d ) = (4.J) 

f p (d) 

It  is  often  easier  to  maximize  a monotonic  function  of  the 
posterior  density,  such  as  the  logaritnm.  Taking  the  log 
of  Eq.  (4.1)  we  have 

Max  Unp(f|d)  = £np  (d  | f ) + £ np  ( f ) -£np  (d ) ] . (4.2) 

f 

Since  the  last  term  8n  p(u)  on  the  right  does  not  depend  on 
t.»  we  neglect  it  in  maximization  with  respect  to  jf.  Thus 
the  MAP  estimate  equation  is  given  by 

3£np (f | d)  32np(d|f)  3£np(f) 

3T  = 37  + 3T  = - (4,3) 


From  Eq.  (2.31),  we  have 


£np  (d  | f ) = 2* 


£n 


(\f.)di/ae‘>fi 


— a~ — 

a(^)l 


(4.4) 


= 23  £n(Af . )-Af  ,-£na-£n[  (— ) ! ] l . 

. ( a li  nt  | 


From  Eq.  (2.10),  we  have 


£np ( f ) = £nKf-j(f-f )TR"1 (f-F)  . 


(4.5) 


\ 


~r’i;  v . 
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Dif ferentiating  the  above  two  equations  individually  with 
respect  to  l we  get 


np (d| f ) 


9f np(f ) 


Substituting  Eqs.  (4.6)  and  (4.7)  into  Eq.  (4.3) 


Taking  the  transpose  on  both  sides  of  the  above  equation 
and  assuming  that  the  covariance  FfX  is  a symmetric  matrix, 
(i.e.  Rf1=(Rf1)T)  we  then  get 


From  Eo.  (4.9),  we  know  if  the  norm  of  R 
1 1 R f 1 1 is  very  large,  then 


denoted 


-MAP  ' - -ML 

/v 

where  d is  the  observation  data  vector  and  f 


maximum  likelihood  estimate  vector.  In  the  blurred  image 


case,  the  J^L  is  the  inverse  solution  instead  of  the 
observation  data. 


On  the  other  hand,  if  the  norm  of  is  very  small 


then 


f ~ f 
—MAP  ~ - 


where  T is  the  a priori  nonstationary  mean  of  the  image. 


Therefore,  solving  these  equations  for  XN1Ap  tries  to 


move  the  solution  of  f from  the  a priori  nonstationary  mean 


t to  a maximum  likelihood  estimate  f,„T  . Here  Rc  is  a 
— —ML  f 


measure  of  our  confidence  in  the  nonstationary  mean  f ano 


maximum  likelihooa  estimate  fWT  as  a solution  to  the 

—ML 


restoration  problem.  Equation  (4.9)  appears  very  simple, 
but  the  complexity  of  the  estimate  implementation  depends 


heavily  on  the  structure  of  the  R^  covariance  matrix.  Thus 


we  will  aiscuss  in  the  two  following  sections  methods  of 
implementing  Eq.  (4.*).  One  assumes  R^  is  an  identity 


matrix,  anci  the  other  assumes  that  Rf  is  a Markovian 


matr  ix . 


4 . 3 MAP  Estimate  Implementatin  with  an  A Priori  Image 
Covar lance  Matrix  - Identity  Covar  iance  Matr ix 


FoR  simplicity,  we  assume  R^=0^I , thus  each  pixel  of 
the  picture  is  uncorrelated.  in  real  picture  data,  each 
pixel  is  highly  correlated  with  its  neighbors  | 4 — 1 1 1 , and 
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-T  v 


■ nr. 


I 


this  assumption  is  treated  in  the  next  section.  From 

2 

Eq.  (4.9)  and  ^^=0^1  we  then  have 


■^-a' 

fl"?l 

r ' 

o 

i 

d 

«f2-* 

1 

- -j 

f2"?2 

= 

0 

dN  . 

_ a f N 

fN_?N 

0 

From  Eq.  (4.10),  we  see  that  the  MAP  estimate  becomes  a 
very  simple  point  process  instead  of  a vector  process 
because  equations  are  decoupled.  Hence  we  can  get  a closed 
form  solution 

(7. - Ao^ ) +yf(f-  -Ao ?)+4Aa?d . 

f = — i £ i - — (4.11) 

1 2 

where  the  positive  root  is  taken  because  intensity  is 
always  non-negative. 

4.3.1  Implementation  and  Exper imental  Results 

The  observation  data  are  photon  counts  with  some 
amplification  gain  a.  The  photon  count  is  simulated  from 
an  original  picture  (256x256)  through  a Poisson  random 
noise  generator  wich  is  described  in  detail  in  Appendix  A. 

In  order  to  implement  Eq.  (4. 11), we  must  estimate  the 

, 2 — 
variance  of  and  nonstationary  mean  J[  from  the  available 

aata.  Hunt  [4-131  has  shown  that  an  estimate  of  the 

nonstationary  ensemble  mean  of  a context-dependent  ensemble 


can  be  made  by  blurring  the  given  ensemble  member  with  a 
point-spread  function  as  expressed  by 

<fK(x,y)>  = f (x,y)  = h (x,y)0f (x,y)  (4.12) 

where  <f  (x,y)>  denotes  an  ensemble  average  at  (xfy)  over 

K 

all  k ensemble  images 

t(x,y)  is  the  nonstationary  mean  of  the  object 
f(x,y)  is  the  object  intensity 
0 denotes  two-aimens ional  convolution 
h(x,y)  is  a point-spread  function  related  to  the 
probabilistic  process  which  generates  the 
ensemble . 

Equation  (4.12)  states  that  an  estimate  of  the  ensemble 
mean  can  be  made  by  convolving  the  original  image  with  a 
point-spread  function  h(x,y).  In  the  experimental  work 
shown  in  this  section,  an  estimate  of  the  nonstationary 
mean  along  each  line  is  made  by  blurring  the  noisy 
measurement  data  with  a point-spread  function  h(x,y)  chosen 
to  be  square  blurring  function.  We  have  chosen  an  11  pixel 
linear  moving  average  window  along  each  line  in  order  to 
obtain  nonstationary  mean  estimate.  Because  of  this  large 
moving  window,  the  noise  is  averaged  out  while  the  low 
frequency  components  of  the  object  remain.  The  image 
variance  is  estimated  by  an  unbiased  estimate  of  the 
population  variance. 


The  restored  pictures  produced  by  the  solution  of 


are  shown  . in 


signal  to  noise  ratio  denoted  by  (SNR) 


it  can  be  seen  that  the  restored  pictures  are  improved 


compared  to  the  noisy  pictures  and  that  the  nonstationary 


mean  carries  the  important  low  frequency  image  information 


4.4  MAP  Estimate  Implementation  With  An  A Priori 


In  this  section 


we  assume 


covariance  matrix  with  correlation  coefficient  P.  The 


Markovian  covariance  matrix  is  very  good  approximation  for 


real  image  signals.  The  Markovian  covariance  matrix  for  a 


one-aimensional  image  model  is 


It  can  be  shown  that  the  inverse 


Restored  imaae  with  a MAP  filter  with 
Rc  = a2j  and' (SNR) =/7 


Original  object  image 
Poisson  noisy  imaqe 
Nonstationary  mean  image 
Restored  image  by  the  MAP  filter 


Figure  4.1b  Restored  image  with  a MAP  filter  with 
R,=a2I  and  (SNR)  =/T0 


Original  object  image 
Poisson  noisy  image 
Nonstationary  mean  image 
Restored  image  by  the  MAP  filter 
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i+p‘ 
r = — 


l-p‘ 


1 

~Z 


0 = 


l+p' 


| p | <1  and  | 0 | < 2' 


Substituting  Eq.  (4.14)  into  Eq.  (4.9),  we  obtain  a set  of 
three  types  of  nonlinear  equations  with  N unknowns 


1 1+p 


d. 

l 


(J^-X)t6r(£i.1-r..1l-r(fi-fi)+6r(£.+1-£.+1)  - 0 


(4.14) 


i = 2 , 3 , . . . ,N-1 


,-N_x)+6r(£N.1-fN.1)--i7(£N-rN)  = 0 . 


Due  to  the  large  dimensionality  and  nonlinear  nature  of  the 
above  system  equations,  ordinary  linear  signal  processing 
matrix  operations  and  Fourier  methods  are  of  no  use.  Thus, 
it  is  impossible  to  directly  solve  those  equations  in  order 
to  obtain  optimal  solutions.  Instead  a suboptimal  method 
with  sectioning  to  reduce  the  dimensionality  of  the 
equations  is  used.  After  trying  several  techniques  for 
numerical  solution  of  the  nonlinear  equations,  we  have 
found  that  a Newton-Raphson  iterative  method  (4-4)  is  the 
best.  This  method  is  very  easily  implemented  and  converges 
rapidly  with  a good  starting  guess.  In  this  application, 
convergence  is  generally  reached  in  about  three  to  four 
iterative  steps.  The  details  of  applying  this  method  to 
the  MAP  estimate  will  be  discussed  in  Appendix  B. 
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The  estimate  of  image  variance  of  and  the 
nonstationary  mean  are  done  by  the  same  technique  as  in  the 
last  section.  The  initial  value  of  f in  the  Newton-Raphson 
method  might  be  assumed  to  be  the  nonstationary  mean  jf  or 
the  raw  observation  data,  but  the  final  estimates  must 
converge  to  the  same  values.  The  convergence  criterion  is 
based  on  the  numerical  closeness  to  the  ideal  image 
(object) . The  convergence  rate  is  also  controlled  by  the 

A 

2 

estimated  variance  af  of  the  image.  An  accurate  estimate 
of  the  local  variance  is  very  important  to  the  convergence 
of  the  algorithm.  The  iterative  Newton-Raphson  r ethod 
employs  the  gradient  of  the  function  to  obtain  the 
increment  value  for  iteration,  thus  it  converges  much  more 
rapidly  than  other  numerical  methods. 

The  boundaries  between  sections  must  be  carefully 
considered  when  the  MAP  equations  are  solved  by 
Newton-Raphson  techniques.  The  overlap  sectioning  method 
may  minimize  the  boundary  effects,  but  more  computing  time 
is  required  because  the  number  of  arithmetic  operations  in 
the  Newton-Raphson  solution  goes  up  roughly  as  the  cube  of 
the  section  size.  Thus,  there  is  a compromise  between 
section  size  and  processing  speed.  The  sectioning 
technique  is  best  explained  with  the  use  of  diagram  shown 
in  Fig.  4.2.  In  implementing  the  one-dimensional  MAP 
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-H  No  K- 

N„/2  Nv  N„/2 


(1)  Nq  measured  image  data  points  are  processed  in  each 
section  with  an  overlap  of  Nq  data  points  between 
sections . 

(2)  Each  section  keeps  Nv  valid  processed  data  points  and 
discards  ND/2  erroneous  processed  data  points  at  the 
two  ends. 

(3)  N_ , Nq,  Nv  and  NQ  are  the  section  size,  overlap  size, 
valid  processed  data  size  and  discarded  data  size 
respectively. 


Figure  4.2  One-dimensional  sectioning 
method  diagram 


filter  we  have  found  that  a section  size  of  32  pixels  and 
overlap  of  16  pixels  in  each  section  works  well.  This 
choice  reduces  the  processing  time  and  also  minimizes  the 
boundary  effects.  These  assumptions  are  based  on  the  fact 
that  pixels  separated  by  16  unit  sampling  distances  are 
essentially  uncorrelated  even  if  p=0.95.  The  restored 
pictures  with  the  MAP  filter  using  different  overlap  sizes 
in  each  section  are  shown  in  Fig.  4.3  and  Fig.  4.4. 
Figure  4.3  is  the  restored  picture  produced  by  an  MAP 
filter  with  b pixels  overlap  in  each  section.  Figure  4.4 
is  the  MAP  restored  picture  with  lb  pixels  overlap  in  each 
section.  From  Fig.  4.3  and  Fig.  4.4,  we  can  see  some 
sectioning  boundary  effects  although  the  restored  picture 
of  Fig.  4.4  has  reduced  edge  effects  compared  to  Fig.  4.3. 
The  sectioning  edge  effects  of  Fig.  4.4  are  almost 
invisible.  Thus,  an  overlap  of  16  pixels  in  each  section 
is  a good  practical  choice  for  minimizing  the  sectioning 
boundary  effects.  Therefore,  all  the  following 
one-dimensional  processed  pictures  have  used  the  overlap 
sectioning  method  with  a section  size  of  32  pixels  and  an 
overlap  of  16  pixels  in  each  section. 

We  have  also  found  that  the  epu  processing  time  of  the 
MAP  filter  with  P=U.95  is  lesser  than  that  of  the  MAP 
filter  with  o-u.U  because  of  the  a priori  knowledge  of 
pixel  correlation.  The  restored  pictures  with  MAP  filters 
are  shown  in  Fig.  4.5  for  different  (SNR)rms  with  p =0  and 


Restored  imaqe  with  a MAP  filter  8 pixels 
overlap  in  each  section 


Original  object  image 
Poisson  noisy  image 
Nonstationary  mean  image 
Restored  image  by  the  MAP  filter 


Restored  image  with  a MAP  filter  16  pixels 
overlap  in  each  section 


Original  object  image 
Poisson  noisy  image 
Nonstationary  mean  image 
Restored  image  by  the  MAP  filter 


Figure  4.5a  Restored  image  with  a MAP  filter  at 

(SNR)  =/5  and  0=0 
rms 

(A)  Original  object  image 

(B)  Poisson  noisy  image 

(C)  Nonstationary  mean  image 

(D)  Restored  image  by  the  MAP  filter 


Restored  image  with  a MAP  filter  at 

(SNR)  =/5  and  p=0.95 
rms 

(A)  Original  object  image 

(B)  Poisson  noisy  image 

(C)  Nonstationary  mean  image 

(D)  Restored  image  by  the  MAP  filter 


Restored  image  with  a MAP  filter  at 

(SNR)  =/T0  and  o=0.95 
rms 

(A)  Original  object  image 

(B)  Poisson  noisy  image 

(C)  Nonstationary  mean  image 

(D)  Restored  image  by  the  MAP  filter 


Restored  imaqe  with  a MAP  filter  at 

(SNR)  =/20  and  o=0 
rms 

(A)  Original  object  image 

(B)  Poisson  noisy  imaqe 

(C)  Nonstationary  mean  image 

(D)  Restored  image  by  the  MAP  filter 
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Restored  image  with  a MAP  filter  at 
(SNR)  =/20  and  p=0.95 


Original  object  image 
Poisson  noisy  image 
Nonstationary  mean  image 
Restored  image  by  the  MAP  filter 
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p =u . 95 . The  restored  images  are  optimal  solutions  between 
the  maximum  likelihood  (ML)  solution  and  the  a pr ior i 
nonstationary  mean  solution.  The  MAP  solution  smooths  out 
the  Poisson  noise  degradations  and  has  also  extracted  some 
higher  frequency  information  from  the  noisy  images.  The 
images  restored  with  the  MAP  filter  have  more  detail 
information  particularly  at  higher  (SNR)rms . 

4 . b Recursive  MAP  Estimate  for  the  Poisson  Noise  Model 

From  Eq.  (4.14)  and  simple  algebraic  manipulations,  we 

get 


— l x ^ 1 

f,  = f9+ ^-(f.-f. )-— (^-l) 

* 8(1+0  ) Br  1 


(4.15a) 


£i+l  ■ fitl+F(£i'ri)-IF(r7‘1,'<fi-l'£i-l)'  «-15b> 


i = 3,4, .. . , N-l 


i 


A 


I 

I 


i 


recursive  relation 


K+l  = fi+kg(fi,di,f)  (4.16) 

where  fi  is  the  previous  estimate,  di  is  the  observation 

A 

data,  k is  a residual  gain,  F is  the  estimate  of  the 

/\  /v 

nonstationary  mean  vector  and  g is  a function  of  f^,  d^,  F. 
Because  images  with  Poisson  noise  have  a very  low  (SNR) rms , 
the  nonstationary  mean  can  only  be  estimated  approximately. 
Also,  the  estimation  error  is  propagated  through  all  the 

A 

estimates  f^.  Thus,  Eq.  (4.16)  is  very  unstable  and  it  is 
impossible  to  obtain  an  accurate  recursive  solution.  Some 
simulations  have  been  performed  with  different  estimates  of 

A A 

f^  ana  f^  for  different  (SNR)rms.  All  experimental  results 
quickly  blow  up,  obliterating  all  the  image  information. 


4.7  A Local  Adapt ive  Processing  Filter 

As  discussed  in  Ch.  2,  image  signals  are  a 

nonstationary  random  field  whose  statistical  properties 
vary  in  a local  region  of  the  image.  Hence,  a local 
adaptive  processing  filter  should  have  many  advantages 
compared  to  global  processing  filters  which  are  defined 
over  the  entire  image  field.  Global  processing  filters 
generally  average  over  detail  information  in  local  regions 
of  the  images.  The  local  adaptive  filter  should  be 

particularly  useful  when  the  image  noise  is 
signal-dependent  as  with  Poisson  noise.  The  sectioned  MAP 
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filter  as  described  can  be  made  to  operate  as  a local 


adaptive  processor.  Because  the  MAP  filter  contains  an  ML 
(maximum  likelihood)  term  and  an  a pr ior i term 
(nonstationary  mean)  , the  filter  can  be  implemented  by 
adaptively  weighting  the  terms  as  a function  of  local 
properties  such  as  the  first  and  second  moments  of  the 
image  or  the  nonlinearity  of  the  human  eye.  The  local 
adaptive  MAP  filter  can  be  expressed  as 

* (ML  term)  + (1-Q^)  * (a  priori  term)  = 0 


where  is  the  weight  of  the  ith  section  and  L is  the 
total  number  of  sections  in  the  entire  image.  The  local 
adaptive  filter  can  also  be  extended  to  space-variant 
degradations.  For  simplicity  in  the  experimental  results 
that  follow,  we  have  simulated  the  global  adaptive 
processing  filter  with  Qj=Qj=Q  for  all  sections.  The 
equations  to  be  solved  are: 


Q*( -r-  1)-(1-Q)*< T-^~)  (f1-f1)  + (l-Q)*X{f2-f2)  (4.18a) 

1 (1+P  )A  _ o 


1 

+ I1-Q >*T  <fi+l-W  - 0 


Q*(/-l)  + (l-Q)  *X(fi-l'?i+l)-(1"Q)* 


V 


(fN'V  " ° 


(1+P2) 


(4.18b) 


(4.18c) 


A Newton-Raphson  iterative  technique  with  sectioning  was 


r 


used  to  solve  Eq.  (4.18).  All  the  Newton-Raphson  iterative 
techniques  are  described  in  Appendix  B.  The  only  factor 
that  changes  are  some  constant  coefficients  in  the 
equations. 

The  resulting  pictures  of  this  global  adaptive  filter 
are  shown  in  Figs.  4.6  through  4.6  for  different  (SNR)rms 
and  different  weights  Q=0.3,  0.5,  and  0.8.  The  Q=0 . 5 
weights  both  terms  equally.  The  region  of  higher  (SNR) rms 
can  be  given  a higher  weight  and  then  can  extract  higher 
frequency  components  from  the  ML  term,  while  the  region  of 
lower  (SNR) rms  can  be  weighted  more  towards  the  a pr ior i 
term.  These  facts  can  be  seen  in  Figs.  4.6  through  4.8 
where  we  can  see  clearly  that  detail  information  is  more 
visible  and  noise  suppression  is  decreased  as  the  ML  term 
is  increased.  However,  if  Q is  too  large,  then  the 
resulting  image  is  the  same  as  the  unprocessed  data. 

4.8  Conclusions 


I 

l 


1 


4 


From  this  chapter,  we  have  found  that  the  estimated 
nonstationary  mean  carries  most  of  the  gross  information  in 
MAP  estimation,  and  that  the  covariance  matrix  carries 
important  detail  information  of  an  image.  Also,  the 
variance  affects  the  convergence  rate  of  the  algorithm. 
The  variance  can  act  as  a weighting  factor  in  sectioned 
suboptimal  MAP  estimation  because  this  method  is  very 
dependent  on  the  local  nonstationary  variance.  The 
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Poisson  noisy  image 
Restored  image  with  w 
Restored  image  with  w 
Restored  image  with  w 


' 


Restored  images  with  an  adaptive  MAP 

filter  at  (SNR)  =/20 
rms 

(A)  Poisson  noisy  image 

(B)  Restored  image  with  w.=0.3 

(C)  Restored  image  with  w.=0.5 

(D)  Restored  image  with  w.=0.8 


sectioning  method  with  a Newton-Raphson  solution  has  been 


effective  at  coping  with  large  dimension  nonlinear  MAP 


estimation  equations  and  has  produced  good  results 


CHAPTER  5 


RESTORATION  OF  BLURRED  IMAGE  WITH  POISSON  NOISE 

5. 1 Introduction 

This  chapter  we  extend  the  previous  results  to  the 
more  general  imaging  model  including  blurring  degradations 
and  Poisson  noise  degradations  as  discussed  in  section  2.4. 
In  many  practical  situations  of  interest,  the  detected 
image  data  arises  from  a linearly  blurred  image  of  the 
object.  Although  the  blurring  may  arise  from  many 
different  sources,  we  simply  lump  them  together  as  a 
blurring  matrix  H.  This  system  and  its  block  diagram  are 
shown  in  Fig.  2.5  and  Fig.  2.6  respectively.  In 

section  5.2,  we  formulate  the  MAP  estimation  equations  and 
its  solution.  In  section  5.3,  we  briefly  review  the 
sampled  infinite  area  superposition  operator  model  for 
image  blur.  In  section  5.4,  we  discuss  sectioning  methods 
and  in  section  5.5,  we  discuss  implementation  of  an  MAP 
filter  with  one-dimensional  blurring  and  its  experimental 
results.  In  section  5.6,  we  discuss  implementation  of  the 
MAP  filter  with  two-dimensional  blurring  and  present  some 
experimental  results.  In  the  final  section,  we  present 
some  conclusions  of  this  chapter. 


■»nf 


60 


5 . 2 MAP  Derivation  with  A Blurring  Matrix  H 

From  Bayes'  Law  we  have  the  posterior  density  as 
follows : 

p (d  | b)  p (b) 

Plbldl ^ (5.1) 

where  b=Hf  and  b.=  E 

1 j 13  3 

From  Fig.  2.6,  we  know  that  cj  is  the  Poisson  noise 
degraded  version  of  (Hf) . As  before  we  have  the 
conditional  density 

g . -Xb . 

(Xb.)  Xe  1 

Pfqjb  ) = (5.2) 

yi* 

for  the  measured  counts  g.  as  a function  of  incident 

1 

intensity  b^.  We  also  assume  that  counters  i and  j are 
independent  for  a given  b.  Hence 

p(2|b)  = p (gx | b) p (g2 | b) p (gN | b)  (5.3) 

Because  each  g^  only  depends  on  its  corresponding  b 

p(2|b)  = p(gx |b1)p(g2 |b2) . . .p(gN|bN) 

g -Xb. 

(Xb.)  xe  1 (5.4) 

= n 

ana  setting  d ,*ag  . 

11 


we  have 


b = H f 


(5.10) 


Substituting  Eqs.  (5.9)  and  (5.10)  into  Eq.  (5.8),  we  get 


p (Hf ) = Kfeexp{-j(f-f )TRf1 (f-f ) } 


(5.11) 


5.2.2  Estimation  Equations 


From  chapter  4,  we  have  the  estimation  equations  as 


f ol lows : 


£np(d|b)+A-£np(b)  = 0T, 


(5.12) 


where 


b = Hf  and  b.  = 53  H. .f. 
- - i j ^ 1 


(5.13) 


From  Eq.  (5.5),  we  have 


Hnp(d|b)  = ^ 


d . /a  - Ab . 

l l 


(Ab.)  1 e 

<T 

a(-i)  ! 


(5.14) 


Ed,  d . 

{ — £n(Ab.)-Ab.-£n«-en[  ( — 


{-^■dn(Abi)-Abi-Hnoi-tnl  (~)  !]  } (5.15) 


8,np  (d  | b)  = 


Ld  . H . 


(5.16) 


V d. 

i A{Bi-HiK) 


A:  t 


(5.17) 


4?*np(d|  Hf)  = X (gi 

L i i 


i i 


(5.18) 


and  from  before 


?x(Ei_1)HiN 
x 1 


^p(Hf)  = - (f-f ) Tr£  1 


(5.19) 


Substituting  Eqs . (5.1b)  and  (5.19)  into  Eq.  (5.13),  we  get 

d, 

-1  )U 

i2 


I ^ X *b^_1*  Hil ' ^ X (57-1*  Hi2  ' • • * ' ^ X ^-1)  H1n1 
L i i ii  ii  J 


+ [-(f-r)'I’Rf1]  = 0T 


(5.20) 


Taking  the  transpose  of  Eq.  (5.20)  and  assuming  Rf=R* 


then,  we  get 


Ed . 
x X(b7'1)Hil 


Ed  . 

X (Ei_1)Hi2 
1 1 


Ld  . 

A(E7“1)HiN 


.-1 


■ Rf  (f-f)  = o (5.21) 


We  expand  the  first  term  of  the  left  side  of  Eq.  (5.21)  to 


obtain 


d,  ^ 

^-1)Hn+V1)H21+* 


N 


. + (b  !)Hn1 

N 


d,  ^2 

^-1)H12+(b^-1)H22+- 


N 


,+ (bN~1)HN2 


d ' d2  dN  ' 

L (b^~1)HlN+(b^  1)H2N+- ' •+ (bN  1)HNN 


1 

'Hll 

H21 

H31  •“  HN1 

H12 

H22 

H32  HN2 

= X 

• 

• 

• 

• 

• 

• 

-hin 

HN2 

HN3  * • ‘ HNN 

H11 

H12 

....  h1n 

= X 

H2  i 

H22 

....  h2n 

_hni 

HN2 

* • • • hnn 

= XHT(q-i) 

i 

Here  the  H 

matr ix 

is  not 

necessarily 

dl 

r r-i-  1 
bl 


- 1 


N 


- 1 


T r 


£-i 

bl 


d2 

y~~  ~ 1 
b2 


?H-1 

bN 


(5. 22) 


(5.23) 


(5.24) 


(5.25) 


depends  on  the  model  of  the  blurting  degradation  where 

i T 


3.  _ 1 3 1 »*!  2 ' * ’ ’ '^N  ^ 


and 


• y 


q 

1 .BJ  ZHijfj 


For  the  time-being,  we  assume  that  H is  a square  matrix 


H11 

H12 

. . . 

1 

2 

rH 

a: 

H2 1 

H22 

• • • 

H2N 

(5.26) 

hni 

HN2 

• • • 

hnn- 

Therefore,  the  MAP  estimate  equations  reduce  to 


XHT(q-l)-Rf1 (f-f)  = O 


(5.27) 


This  equation  is  the  key  equation  of  MAP  estimation  with  a 
blurring  matrix. 

The  complexity  of  solving  this  equation  is  determined 

by  the  structure  of  the  blurring  matrix  H as  well  as  the 

covariance  matrix  of  image  denoted  by  Rf . If  the  matrix  Rf 

is  a separable  matrix  it  can  be  expressed  as  R =R  ©R 

L L R I C 

where  0 denotes  the  direct  product  {5-14]  and  R and  R 

f R f c 

are  NxN  matrices  of  the  form 

r . 2 N-l-i 

1 p p . . . . p 

. N-2 

p 1 p . . . . p 

R£„  = . (5.28) 


N-l  N-2  N-3 

P P P 


The  Markovian  covariance  matrix  R, 


is  an  accurate 


approximation  to  the  statistics  of  many  images  [5-2],  The 
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H matrix  is  an  N xN  blurring  matrix  which  may  be 

nonseparable  or  separable.  If  H is  a separable  matrix, 

then  H=H  «H  where  H or  H is  an  NxN  matrix.  The  detail 
R c R c 


of  modeling  continuous  superposition  integrals  for  blurring 
by  a discrete  operator  H is  discussed  in  the  next  section. 

Solving  the  MAP  estimation  Eq.  (5.27)  is  heavily 
dependent  on  the  separability  of  the  H matrix  and  Rf 
matrix.  If  we  assume  that  both  are  separable,  then  the  MAP 
estimate  Eq.  (5.27)  can  be  obtained  by  separate  row  and 
column  processing.  Thus,  the  computation  time  of  a 
separable  two-dimensional  MAP  filter  is  twice  that  of  a 
one-dimensional  MAP  filter.  However,  the  computation  time 
of  a nonseparable  two-dimensional  MAP  filter  is 
approximately  the  square  of  the  time  for  a one-dimensional 
MAP  filter.  Thus,  computing  H (3-I)  takes  approximately 
2N 4 operations  (each  operation  includes  one  mul t ipl icat ion 
and  one  addition)  for  the  nonseparable  case  and  4N ^ 
operations  for  the  separable  case  where  N is  the  picture 
size,  e.g.  N=256.  The  tremendous  amount  of  computing 
needed  for  the  MAP  estimate  makes  solution  impossible  even 
with  the  separable  case.  Thus,  we  adopt  a suboptimal 
solution  involving  sectioning  with  a Newton-Raphson 
solution  technique.  This  method  has  been  developed  in 
Chapter  4 for  solving  nonlinear  MAP  estimate  equations  of 


larger  dimensionality.  Using  the  sectioning  method  [5-4], 

I 


the  MAP  estimate  Eq.  (5.27)  becomes 


X«T  (q-l)(m)-Rf1  (f-f ) (m)  =0  (5.29) 


where  the  superscript  m denotes  the  mth  section  of  the 

image.  If  each  section  size  is  N 2,  then  3 and 

t are  N2  xl  vectors  and  H ^ is  an  N 2 xN  2 matrix.  For 
— s s s 

clarity,  we  omit  the  superscript  (m)  in  the  discussion 
which  follows.  The  H matrix  is  constructed  from  the  known 
point  spread  function  (PSF)  h(x,y).  The  PSF  h(x,y)  is 
obtained  either  from  a pr ior i knowledge  or  a posterior i 
knowledge.  The  former  assumes  that  the  PSF  h(x,y)  is  known 
a E£i£li-  The  latter  assumes  that  the  PSF  h(x,y)  is  not 
available  a pr ior i . It  must  be  estimated  from  noisy 
observation  data  using  techniques  such  as  blind 
deconvolution.  References  15-5,5-6,5-7,5-8)  describe 


several  methods  for  estimating  the  amplitude  response  of 
h(x,y),  and  the  Knox-Thompson  method  can  be  used  to 
estimate  the  phase  response  of  h(x,y)  [5-9,5-15].  These 
methods  determine  the  degradation  parameters  by  a 
posteriori  methods  (5-3).  This  thesis  assumes  that  PSF 
h(x,y)  of  the  system  is  given  as  a pr ior i information.  We 
must  know  the  detailed  structure  of  the  H blurring  matrix 
in  order  to  implement  the  MAP  estimate  Eq.  (5.29).  Hence, 
we  now  discuss  how  to  construct  the  H matrix  from  a given 
PSF  h(x,y)  of  the  system. 


4 
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5 . 3 Construction  o£  the  Blurring  Matrix  H 

All  blurring  degradation  effects  are  lumped  together 
as  a two-dimensional  point  spread  function  h(x,y;£,n)  as 
shown  in  Fig.  5.1.  From  linear  system  theory,  we  can 
describe  Fig.  5.1  by  the  superposition  integral 


G (x,y) 


F (a,  B)h (x,y;a, B)dotdB  (5.30) 


If  this  linear  system  is  spatially  invariant,  then 

h(x,y;a,B)  = h(x-a;y-B).  (5.31) 

In  order  to  discretize  the  continuous  Eq.  (5.30),  we  must 
sample  G(x,y)  over  a grid  at  spacing  (Ax, Ay)  satisfying  the 
Nyquist  criterion  (5-2J.  For  notational  simplicity,  the 
continuous  object  function  F(ot,B)  and  the  continuous  PSF 
h(x,y;£,n)  will  also  be  assumed  to  be  sampled  over  the  same 
grid  spacing.  Thus,  Eq.  (5.30)  can  be  expressed  as  a 
double  summation  over  infinite  limits  by  invoking  the 
sampling  theorem  and  a quadrature  formula.  In  order  to 
express  the  infinite  area  superpositon  operator  as  a 
lexicographic  ordering  for  vector  processing,  it  is 
necessary  to  truncate  the  PSF  to  some  spatial  limit,  say 
( LAx , LAy ) and  restrict  G(x,y)  to  an  area  (MAx,MAy).  Then 
the  truncated  superposition  operator  is 
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L+m^-1  L+n^-l 

G(m1,m2)  = 53  53  F(n^,n2)p(mj-n^+L,  (5.32) 

n^-m^  r\2”^^2 

m 2_ ni ^ / ^2  ^ 

where  for  simplicity,  the  grid  spacing  Ax, Ay  is  dropped  and 
p,  assumed  to  be  zero  outside  its  range  of  indices, 
represents  the  sampled  point-spread  function  which 
incorporates  the  quadrature  integration.  The  detailed 
derivation  of  Eq.  (5.32)  is  in  [5-2].  In  order  to  prevent 
serious  approximate  error  at  the  boundaries  of  G,  we  must 
choose  N such  that 


N > M+L-l  (5.33) 

where  N is  the  size  of  the  convolved  image  signal  G 

M is  the  size  of  the  object  image  signal  F 

L is  the  size  of  PSf  P 

These  boundary  problems  are  important  and  closely  related 
to  the  sectioning  method  used.  In  the  next  section  we 
discuss  the  use  of  an  overlap-save  sectioning  method  to 
minimize  the  boundary  edge  effects. 

If  the  arrays  F and  G are  represented  as  lexicographic 
ordered  vectors  by  f^  and  3 respectively,  then  the 
superposition  operator  can  be  written  as 

3 = H f (5.34) 

2 2 

where  f^  and  3 are  M xl  and  N xl  vectors  respectively.  H is 
the  xN^  matrix 
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H = 


-11^  -12n 

V \ 


— 1L 

V 


-22 


v —2 , L— 2 , L+ 1 

N \ \ 
H. ,'...  H 


--M,N-L+1 

where  H—  are  MxN  matrices  with  entries 


--M , N 


(5.35) 


H (m. ,n. ) = p (m. -n. +L,nu-n0+L) 

m2  ^2  11  11  22 


(5.36) 


for 


i 


1 <_  <_  M, 


1 — m2  — M 


ml  1 ni  — m2  — n2  — L+m2-^ 


(5.37) 


and  P is  a discretized  truncated  point  spread  function 
h(x,y).  if  the  PSF  is  spatially  invariant  (SIPSF)  , then 


H = H .. 

— m0,n0  — m2+l,n0+l 


(5.38) 


2 9 2 2 9 2 

When  the  PSF  is  spatially  invariant  and  orthogonally 
separable 


H " Hc  • HR 


(5.39) 


and  the  two-dimensional  convolution  operation  becomes 


G = HcFHr 


(5.40) 


where®  denotes  direct  product,  G and  F are  the  image  and 
object  arrays,  respectively,  and  H_  and  H are  MxN  matrices 

K C 

of  the  form 
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hr  - 


pr(l> 


Pr(L-1) 


Pr(1) 


VL> PR(1)-J 


(5.41) 


Another  discrete  operator  models  the  blurring  degradation 
by  a finite  area  superposition  operator  expressed  as  a D 
matrix.  The  form  of  this  model  is 


Q (m^ ,m2) 


F (n  j , n2 ) p (mi_n1  + l ,m2-n2  + l ) 


(5.42) 


where  M=N+L-1 . Hence  the  processed  array  C is  of  larger 
dimension  than  object  data  array  F.  Its  vector  form  is 


g = D f 


(5.43) 


2 2 

where  D is  an  M xN  matrix  of  the  form 

* °1  1 

' ( 

-2,1  -2,2^ 


D = 


-L,1>^l-1,2  *” 


-L,2. 


N 


(5.44) 


where  D — is  an  MxN  matrix. 


All  the  special  cases  of  the  D matrix  are  the  same  as 
those  of  the  H matrix  except  for  the  matrix  structure  which 
has  the  form  of  Eq.  (5.44).  The  difference  between  the  H 
operator  matrix  and  the  D matrix  for  modeling  a continuous 
superpositon  are  that  the  processed  array  for  the  finite 
area  H computation  is  equivalent  to  the  processed  array  for 


..  .■ 
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the  finite  urea  H computation,  surrounded  by  a boundary  of 
L-l  superfluous  elements.  Conversely,  if  the  processed 
array  size  is  the  same  for  two  superposition  operators,  the 
L-l  boundary  elements  for  the  array  obtained  by  the  D 
operator  will  be  in  error.  The  resulting  implementation  of 
the  MAP  filter  using  the  D operator  for  the  blurring  matrix 
is  shown  in  Fig.  5.2.  Figure  5.3  shows  the  results  for  the 
H operator.  We  can  see  the  checkerboard  noise  appearing  in 
the  restored  image  of  Fig.  5.2.  The  checkerboard  noise 
results  from  using  the  D operator  giving  rise  to  the 
erroneous  data  surrounding  each  section.  Thus,  the 
two-dimensional  MAP  filter  is  implemented  with  an  H 
operator  matrix  for  the  nonseparable  blurring  cases  rather 
than  the  D operator.  We  can  either  add  zero  elements  to 
the  rectangular  H matrix  or  truncate  the  boundary  elements 
which  wrap  around,  making  it  a square  matrix  as  discussed 
in  previous  sections.  However,  the  discrete  blurring 
matrix  H should  be  a rectangular  matrix  physically.  Since 
m finite  data  points  are  convolved  with  n finite  data 
points,  the  resulting  processed  data  has  m+n-1  points. 


5.4  Sectioning  Method 


Due  to  the  large  dimensionality  and  nonlinearity  of 
the  MAP  estimate  equations,  a sectioning  method  15-4,5-10] 
is  used  with  the  Newton-Raphson  solution  to  obtain  a 
suboptimal  solution.  There  are  two  sectioning  methods: 
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Figure  5.2  Restored  image  with  the  MAP  filter  with  a 
D operator  for  modeling  the  blurring 
degradation 

(A)  Original  object  image 

(B)  Degraded  image 

(C)  Nonstationary  mean  image 

(D)  Restored  image  by  the  MAP  filter 


Restored  imaqe  with  the  MAP  filter  with  an 
H operator  for  modeling  the  blurring 
degradation 

(A)  Original  object  image 

(B)  Degraded  image 

(C)  Nonstationary  mean  image 

(D)  Restored  image  by  the  MAP  filter 
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overlap-add  sectioning  and  overlap-save  sectioning.  The 
sectioning  method  chosen  will  give  rise  to  various  boundary 
edge  effects.  Hence,  we  will  investigate  which  method  is 
the  best  for  MAP  estimation.  From  Eq.  (5.29),  we  have 


outputs  of  the  mth  and  m+lth  section  are  added  together  in 
the  region  of  overlap  to  create  the  final  output.  This 
method  is  suitable  only  for  linear  estimation. 
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an  adjacent  section  m+1 , then  the  overlap-add  sectioning 

method  is  not  valid  because  of  the  nonlinear  function  in 

the  MAP  estimate  equation.  Fortunately,  the  overlap-save 

method  remains  valid  for  the  nonlinear  case  and  can  be  used 

in  our  MAP  estimate  equation  because  incorrect  boundary 

points  in  the  overlap  region  are  discarded,  rather  than 

being  corrected  by  addition.  Thus  the  overlap-save  method 

can  reduce  the  boundary  edge  effects  because  it  discards 

the  erroneous  processed  data  of  the  overlapped  region. 

Since  Eq.  (5.29)  contains  two  H operators  and  assumes  that 

the  truncated  point  spread  function  matrix  is  LxL,  the 

amount  of  overlap  required  is  2 (L-l ) x2 (L-l ) . If  we  must 

correct  NxN  points  of  processed  data  at  each  section,  then 

we  must  use  a working  section  of  (N+2 (L-l ) ) x [N+2 (L-l ) ] . 

The  necessary  overlap  area  constitutes  the  major  overhead 

in  the  sectioned  filtering  process.  It  is  clear  that 

smaller  sections  have  a larger  percentage  of  overhead 

computation.  It  is  also  clear  that  the  computation  will  be 

more  inefficient  for  larger  point  spread  functions. 

However,  it  should  be  kept  in  mind  that  the  number  of 

arithmetic  operations  in  the  Newton-Raphson  solution  for 

each  section  is  the  key  computing  burden  of  the  sectioning 

filter.  In  order  to  find  the  updated  incremental  estimate 

vector  in  each  iterative  step  of  the  Newton-Raphson  method, 

we  must  solve  a set  of  linear  system  equations.  The 

2 

dimension  of  the  linear  system  equations  is  Ns  where  Ns  is 
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the  number  of  pixels  in  each  section  and  the  structure  of 
the  gradient  matrix  depends  on  the  structure  of  the 
blurring  matrix  H when  the  covariance  matrix  R^  is  assumed 
first  order  Markovian.  Hence,  the  number  of  arithmetic 
operations  for  solving  linear  system  equations  is  heavily 
dependent  on  the  form  of  the  blurring  matrix  H.  In 
general,  the  number  of  arithmetic  operations  in  solving  a 
set  of  linear  system  equations  goes  up  roughly  as  the  cube 
of  the  order  of  the  system  when  the  gradient  matrix  is  a 
general  square  matrix.  Thus,  the  smaller  the  section  size, 
the  less  the  computing  time  of  sectioned  filtering  with  the 
Newton-Raphson  solution. 

5.5  Implementation  of  the  MAP  Filter  with  One-dimensional 
Blur r inq  Degradation 

As  stated  earlier,  two  of  the  most  interesting  sources 
of  blur  are  atmospheric  turbulence  and  linear  motion 
[5-11,5-12].  In  this  section,  we  assume  that  the  image  is 
degraded  by  one  dimensional  linear  motion  blur  with  a 
rectangular  point-spread  function.  The  rectangular 
blurring  degradation  is  troublesome  because  its  amplitude 
response  has  a singularity  and  phase  reversals.  The 
blurring  matrix  H is  the  form  of  Eq.  (5.41),  where  pR(L)  is 
the  discrete  truncated  point  spread  function  of  h(x,y).  We 


assume  that  Rf  is  a first  order  Markovian  covariance 
matrix,  and  following  Eq.  (5.29),  we  can  write  the 


equations  to  be  solved  as 
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C2  = 


(1+P2)  (1-P2) 


C1  = ¥ = pC2 

N is  the  number  of  pixels  of  each  section, 

p is  the  correlation  coefficiency  between  pixel, 

2 . 

af  is  the  variance  of  the  object, 

is  the  ijth  element  of  the  H matrix, 

F.  is  the  nonstationary  mean  of  the  ith  pixel  of  the 
section , 

dj  is  the  observation  measurement. 

Using  the  overlap-save  sectioning  method  with  an 
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iterative  Newton-Raphson  solution,  the  MAP  solution  to 


Eq.  (5.47)  is  obtained.  The  convergence  is  very  rapid, 
generally  requiring  about  two  to  three  iterative  si^ps  as 
described  in  detail  in  Appendix  B.  The  discrete  point 
spread  function  of  h(x,y)  is  a rectangular  blurring 
degradation  with  a width  of  5 pixels.  The  nonstationary 
mean  is  estimated  by  a one-dimensional  moving  average  on  11 
pixels  of  observation  data  and  its  variance  is  obtained 
from  all  the  picture  data  by  an  unbiased  estimate.  The 
linear  system  of  equations  for  the  gradient  function  of 
Eq.  (5.47)  which  determines  the  increment  value  for  the 
iterations  is  heavily  dependent  on  the  structure  of  the 
blurring  matrix  H.  When  the  H matrix  is  symmetrical,  the 
computing  time  of  Eq,  (5.47)  with  the  Newton-Raphson 
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Figure  5.4a 


Restored  imaqe  with  a MAP  filter  for 
one-dimensional  blurring  degradation  at 
(SNR) =/5 


rms 


(A)  Original  object  image 

(B)  Degraded  image 

(C)  Nonstationary  mean  image 

(D)  Restored  image  with  MAP  filter 


Restored  image  with  a MAP  filter  for 

one-dimensional  blurring  degradation  at 

(SNR)  =/7T5 
rms 

(A)  Original  object  image 

(B)  Degraded  image 

(C)  Nonstationary  mean  image 

(D)  Restored  image  with  MAP  filter 


Figure  5.4c 


Restored  image  with  a MAP  filter  for 
one-dimensional  blurring  degradation 


(SNR) 


rms 


= /T0 


at 


(A)  Original  object  image 

(B)  Degraded  image 

(C)  Nonstationary  mean  image 

(D)  Restored  image  with  MAP  filter 
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For  a local  adaptive  MAP  filter,  the  equation  to  be 
solved  is 

W- (XHT(q-l) ]+(I-W) • (-Rf1 (f-f ) 1 =0,  (5.48) 

where  W=  diag{W^}.  Here,  is  the  weight  of  the  ith 
secton  which  varies  with  the  nonstationary  mean  and  the 
second  moment  of  the  local  properties  of  the  image.  The 
local  adaptive  MAP  filter  also  can  be  used  for  the 
restoration  of  images  degraded  by  spatially  variant  point 
spread  functions.  The  point  spread  function  at  each  photon 
detector  may  not  be  identical  over  the  whole  array.  The 
image  can  be  divided  into  sectioned  images  each  with  its 
own  space  invariant  PSF.  For  simplicity,  we  have  simulated 
a global  adaptive  MAP  sectioning  filter  in  which  each  W is 
fixed.  The  simulation  is  performed  with  p=0.95  for 
different  section  weights  and  different  (GNR)rms.  The 
experimental  results  are  shown  in  Figs.  5.5  through  5.7 
with  different  (SNR) rms • The  W =0.5  gives  equal  weight  to 
the  (ML)  solution  and  the  a priori  solution. 

From  these  experimental  results,  we  observe  that  more 
high  frequencies  are  extracted  if  the  weight  on  the  ML  term 
is  increased.  Also  we  see  that  large  weight  on  the  ML 
solution  results  in  ill-conditioning  of  some  of  the 
solutions  with  the  MAP  estimate.  With  a large  weight  on 
the  ML  part  of  the  solution,  the  MAP  estimate 
asymptotically  approaches  the  ML  estimate.  This  is 
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Restored  images  with  a global  adaptive 
MAP  filter  for  (SNR) =/5 


rms 


Figure  5.6  Restored  images  with  a global  adaptive 
MAP  filter  for  (SNR)  >/77T 


Degraded  image 
Restored  image  with  w 
Restored  image  with  w 
Restored  image  with  w 
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Restored  images  with  a global  adaptive 
MAP  filter  for  (SNR)  =/l0 


consistent  with  fact  that  the  ML  estimate  is  the  inverse 
filter  for  image  restoration  with  blurring  degradations. 
(For  the  non  blurring  case,  the  inverse  solution  of  image 
restoration  just  is  the  observation  data  as  discussed  in 
the  last  chapter).  From  those  observations,  we  can 
conclude  that  there  is  an  optimal  section  weight  over  the 
global  adaptive  filter  as  well  as  the  local  adaptive 
filter. 

5.6  Implementation  of  the  MAP  Filter  with  a 
Two-d imensional  Blur r ing  Degradation 


This 

section 

is  divided  into 

two 

subsections 

to 

separately 

discuss 

the  assumptions 

of 

separabil ity 

and 

non-separability. 

The  blurring  degrad 

atio 

n is  simulated 

by 

a 3x3  pixel  moving  window  blurring  for  different  (SNR)rms . 
The  nonstationary  mean  is  estimated  by  a 7x7  window  moving 
average  over  the  measured  photon  counts.  This  size  of 
moving  wndow  was  found  to  give  a reasonably  good  estimate 
with  a minimum  amount  of  computing. 

b.6.1  Separable  Case 

In  this  section  references  to  separability  means  that 
the  PSF  is  a separable  space- invar iant  function  (SSIPSF)  in 
the  sectioned  MAP  filter  and  that  the  covarance  matrix  Rf 
is  Markovian  and  separable  [5-3].  This  can  be  expressed  in 
vector  notation  as 
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and  $ denotes  the  direct  product. 


From  direct  product  identities  [5-14],  we  have 

"f1  - (RR  ® V'1  ‘ “r1  ® Rc1  (5.5(1) 

Therefore,  the  two-dimensional  Eq.  (5.29)  can  be 
implemented  using  Eq.  (5.47)  as  a row  processor  and  then 
using  Eq.  (5.47)  as  a column  processor.  The  solution  to 
Eq.  (5.47)  uses  the  same  sectioning  method  with  a 
Newton-Raphson  iterative  solution  as  before.  The 
processing  time  of  this  two-dimensional  MAP  filter  is  twice 
that  of  the  one-dimensional  case.  Without  describing  the 
details  of  the  MAP  implementaton , the  experimental  results 
are  shown  in  Fig.  5.8  for  different  (SNR)rms  and  p=0.95. 

From  Fig.  5.8,  we  can  see  that  the  restored  image  with 
a two-dimensional  separable  filter  performs  considerable 
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Restored  images  with  a separable  MAP 
filter  for  (SNR)  rms=,/5 

(A)  Original  object  image 

(B)  Degraded  image 

(C)  Restored  imaqe  by  one-dimensional 
MAP  filter 

(D)  Restored  image  by  two-dimensional 
MAP  filter 


Figure  5.8a 
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Figure  5.8b  Restored  images  with  a separable  MAP 
filter  for  (SNR)  rms=,/7  • 5 

(A)  Original  object  image 

(B)  Degraded  image 

(C)  Restored  image  by  one-dimensional 
MAP  filter 

(D)  Restored  image  by  two-dimensional 
MAP  filter 
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Figure  5.8c  Restored  images  with  a separable  MAP 

filter  for  (SNR)  =/To 
rms 

(A)  Original  object  image 

(B)  Degraded  image 

(C)  Restored  image  by  one-dimensional 
MAP  filter 

(D)  Restored  image  by  two-dimensional 
MAP  filter 
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smoothing  of  the  Poisson  noise.  The  restored  images  show 
some  improvement  over  the  noisy  originals.  Although  the 
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separable  assumption  is  an  accurate  first-order 
approximation  for  the  system,  the  image  field  itself  is  not 
separable.  Thus,  we  implement  a two-dimensional 
non-sepa rable  MAP  sectioned  filter  next. 

5.6.2  Non-separ able  Case 

We  now  assume  that  the  PSF  is  a non-separable 
space-invariant  function  (NSIPSF)  and  that  Rf  is  an 
identity  matrix  to  simplify  the  simulation.  The 
two-dimensional  non-separable  sampled  infinite  area 
superposition  model  is  used  to  reduce  the  two-dimensional 
blurring  degradation  to  a matrix  H.  This  matrix  is  an 
M xN  matrix,  where  M is  the  observed  data  size  and  N is 
the  processed  data  size  for  the  sectioned  MAP  filter. 
Thus,  we  need  to  solve  a linear  system  of  equations  of 
order  N2xN2  (Eq.  (5.29))  in  order  to  find  the  updated 
increment  value  of  the  root  in  each  iterative  step.  Even 
though  the  sections  are  small,  a large  amount  of  computing 
is  needed  for  this  sectioned  MAP  filter.  The  details  of 
computing  time  of  the  sectioned  MAP  filter  are  described  in 
the  last  section.  Since  the  sectioned  MAP  filter  uses  a 
lexicographic  ordered  vector  representation  and  the 
observed  image  uses  a matrix  representation,  some 
conversion  between  them  is  needed.  The  conversion  relation 
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from  matrix  to  vector  is 


D = ( (J-l) *M) +1  (5.51) 

and  from  vector  to  matrix  they  are 

I = Mod (P-1 ,M) +1  (5.52) 

J=  [ (P-l)Mod(P-l,M) J/M+l  (5.53) 

where  P is  the  location  of  the  vector  element  with 

lexicograhic  ordering,  M is  the  size  of  the  section  to  be 

processed,  (J,I)  is  the  (row,  column)  location  of  the  image 

pixel  and  Mod  is  the  modulo  operator.  The  overlap-save 

sectioned  MAP  filter  is  used  to  minimize  the  boundary  edge 

effects  of  sectioning.  The  simulation  is  done  with  section 

size  9x9  and  an  overlap  of  4x4.  Since  the  blurring 

degradation  is  an  unweighted  average  over  3x3  pixels,  from 

the  discussion  of  the  previous  section,  the  wrap  around 

data  is  2 (L-l ) x2 (L -1 ) pixels  which  is  4x4  pixels.  The 

nonstationary  mean  is  estimated  by  a rolling  window  moving 

average  method  which  efficiently  keeps  only  the  data 

required  in  high  speed  memory.  Becuase  of  the  large 

estimated  CPU  time  for  this  filter  on  the  DEC  KL-10  only 

the  right  side  of  the  noisy  picture  is  processed.  Two 

original  images  with  different  (SNR)  were  filtered  and 

rms 

the  experimental  results  are  shown  in  Fig.  5.9.  From 
Fig.  5.9,  we  can  see  that  the  two-dimensional  filter  has 


Figure  5.9a  Restored  image  with  a non- separable  two- 

dimensional  MAP  filter  for  (SNR)  =/7 . 5 

rms 

(A)  Original  object  image 

(B)  Degraded  image 

(C)  Nonstationary  mean  image 

(D)  Restored  image  by  the  MAP  filter 
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Restored  image  with  a non-separable  two 

dimensional  MAP  filter  for  (SNR)  =/To 

rms 

(A)  Original  object  image 

(B)  Degraded  image 

(C)  Nonstationary  mean  image 

(D)  Restored  image  by  the  MAP  filter 


produced  noticeably  better  results  than  the  restored  images 
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of  Fig.  5. b.  However,  the  CPU  time  of  the  general 
non-separable  filter  is  approximately  a factor  of  100 
longer  than  that  for  a separable  filter  with  a 256x256 
image  and  section  sizes  indicated.  The  results  shown  in 
Fig.  5.9  were  done  with  a preliminary  version  of  the 
algorithm  to  show  feasibility  only.  Processing  the  entire 
256x256  image  would  take  approximately  100  minutes.  It  is 
likely  that  considerable  savings  in  computer  time  would 
result  from  a very  carefully  written  algorithm  which  would 
recursively  perform  Newton-Raphson  solutions  between 
windows . 

5 . 7 Conclusions 

We  have  developed  an  MAP  filter  for  a Poisson  noise 
model  with  blurring  degradations.  The  implementation  and 
method  of  solution  for  the  MAP  filter  are  heavily  dependent 
on  the  form  of  the  blurring  degradation  matrix  H and  the 
covariance  matrix  of  the  object  image.  The  overlap-save 
sectioning  method  with  a Newton-Raphson  solution  has  been 
shown  to  be  an  effective  fast  approach  to  the  suboptimal 
MAP  estimate.  The  sampled  infinite  area  superposition 
model  is  used  for  the  blurring  degradation.  Both  the 
one-dimensional  blurring  and  two-dimensional  blurring 
situations  with  different  levels  of  Poisson  noise  were 
simulated. 
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From  these  experimental  results,  we  find  that  the 
nonstationary  mean  carries  the  most  of  the  structured 
background  information  and  that  the  covariance  matrix  leads 
to  a stable  Newton-Raphson  solution  especially  for  higher 
(SNR)rms  image  signals.  it  has  been  found  that  there  is  an 
optimal  weight  in  the  global  adaptive  MAP  filter  which 
produces  the  best  quality  restoration.  Too  much  weight  on 
the  ML  term  solution  will  gives  rise  to  ill-conditioning. 
From  Fig.  5.9,  we  see  that  the  quality  of  the  restored 
image  with  the  MAP  filter  assuming  a non-separable  PSF 
gives  better  results  than  the  separable  MAP  filter. 
However,  the  CPU  time  for  the  non-separable  filter  is  much 
longer  than  for  the  separable  filter.  The  overlap-save 
sectioned  MAP  filter  has  been  shown  to  be  useful  for 
overcoming  problems  of  large  dimensionality  and 
nonlinearity  which  are  inherent  in  MAP  estimation. 


CHAPTER  6 


THE  QUALITY  OF  THE  MAP  ESTIMATE 

6 . 1 Introduction 

In  this  chapter,  we  investigate  the  quality  of  the  MAP 
estimate  for  Poisson  noise.  The  quality  of  the  estimate 

depends  on  the  performance  criteria  chosen.  There  are  two 

types  of  performance  criteria  [6-1,6-21*.  one  depends  on  the 
estimator  structure  and  the  other  depends  on  the 
performance  itself.  The  MAP  estimate  and  ML  estimate 
belong  to  the  former  one  because  the  MAP  estimate  is  the 

mode  of  the  posteror  probability  density  and  the  ML 

estimate  is  the  mode  of  the  a pr ior i probability  density. 
The  Bayes  estimate  belongs  to  latter  because  it  minimizes 
the  risk  of  the  estimate.  The  MMSE  (minimum  mean  square 
error)  estimate  is  a special  case  of  the  Bayes  estimate 
when  the  cost  function  is  proportional  to  the  mean  square 
error.  However,  it  is  customary  to  choose  the  conditional 
or  unconditional  expected  squared  error  as  a universal 
measure  of  quality  of  all  estimates.  Unfortunately,  the 
expectation  operator  leading  to  this  measure  is,  in 
general,  very  complicated  due  to  the  complexity  of  variance 
estimates.  However,  it  is  possible  to  derive  an  expression 

120 


. ~ a- 


for  a lower  bound  on  the  variance  in  terms  of  only  the 
statistical  properties  of  the  observed  signal  and 
estimation  bias.  This  quality  measure  can  be  formed  for 
any  estimator  without  detailed  knowledge  of  its  structure 
as  long  as  the  estimate  is  unbiased.  This  lower  bound  for 
the  estimation  error  variance  is  the  well  known  Cramer-Rao 
lower  bound  (CRLB)  16-3,6-15]. 

There  are  two  measures  that  are  used  together  to 
determine  the  quality  of  an  estimate.  These  are  the 
expectation  of  the  estimate  and  the  variance  of  the 
estimation  error.  The  first  of  these  measures  the  bias 
inherent  in  an  estimate,  and  the  second  is  equivalent  to 
the  mean-squared  error  between  the  estimate  and  the 
original  data.  In  general,  we  try  to  find  unbiased 
estimates  with  small  estimation  error  variance. 

In  section  6.2  and  6.3,  we  discuss  biased  and  unbiased 

A 

estimates  and  show  that  the  f^p  estimate  for  the  Poisson 
noise  model  is  an  unbiased  estimate.  In  sections  6.4,  6.5 
and  6.6,  we  derive  the  Cramer-Rao  lower  bound  of  the 
estimate  for  non-random  scalar  parameters  and  random  vector 
parameters  with  the  Poisson  noise  model.  In  the  final 
section,  we  present  the  conclusions  of  this  chapter. 

6.2  Biased  and  Unbiased  Estimates 


A conditional  unbiased  estimate  is  one  whose  expected 


value  is  equal  to  the  true  value  of  the  quantity  being 
estimated.  An  unconditional  unbiased  estimate  is  one  whose 
expected  value  is  equal  to  the  expected  value  of  the 
quantity  being  estimated.  We  denote  the  estimate  by  a 

A 

random  variable  X which  is  a function  of  the  observations 
Y.  If  X is  a conditional  unbiased  estimate,  then 


EY  [XI  = 


X (Y) P (Y I X)dY  = X. 


(6.1) 


and  if  X is  an  unconditional  unbiased  estimate,  then 


Eym  = 


X(Y)P(Y)dY  = E (X)  = X. 


(6.2) 


On  the  other  hand,  biased  estimates  do  not  possess  this 
desirable  feature;  their  expected  values  contain  an 
additional  function  B(X)  of  the  parameter  to  be  estimated. 
Accordingly,  for  biased  estimates  we  have 


Ey  [X]  = X+B  (X)  , 


(6.3) 


or 


Ey [X]  = X+bTxT.  (6,4) 

for  the  conditional  biased  estimate  and  the  unconditional 
biased  estimate,  respectively  (6-3,  6-4,  6-5], 

A 

6.3  f map  — an  Unconditional  Unbiased  Estimate  Vector 
From  previous  chapters,  we  have 

ImAP  - rtXRfHT(a-l)  (6.5) 
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where  l^p  is  the  N xl  estimate  vector 
- 2 

f_  is  the  N xl  nonstationary  mean  vector 

Rj  is  the  covariance  matrix  of  the  image 
2 2 

and  H is  the  N xN  discrete  blurring  matrix. 

Taking  the  expectation  on  both  sides  of  Eq.  (6.5),  we  have 


E[fMAP]  = f+Et^RfHX(q-l)]. 


(6.6) 


Since 


d. 

l 

qi  = b~r 

i 


qi  = xbi' 


(6.7) 


then 


E [q±] 


ag^ 

= E,  { E (-r-i)  } 

bi  gilbi  bi 


*“Eb.{Eq.|b.  'sr11 

l I l l 


(6.8) 


where  E |b  denotes  conditional  expectation  over  g^  for  a 


given  b. . E,  denotes  expectation  over  b- , hence 
i 

Xb. 

Etq^  = otEb  [-g-i]  = aX  = 1 


(6.9) 


Thus 


E[q)  = 1 


(6.10) 


ana  substituting  Eq.  (6.10)  into  Eq.  (6.6),  we  get 


E[W  = f- 


(6.11) 


vJ4;  s 
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is  an  unconditional  unbiased  estimate 


r 


I 

I. 


4 

r 


Therefore,  f^p 
vector . 


6.4  Cramer-Rao  Lower  Bound  (CRLB) 

For  notational  and  mathematical  simplicity,  we  first 
focus  our  attention  on  the  CRLB  for  the  case  of  non-random 
scalar  parameters,  an  example  is  sample  mean  and  sample 
variance  of  the  amplitude,  phase,  and  Doppler  frequency  of 
the  estimation  of  the  radar  signals  assuming  these 
parameters  unknown  but  not  random  variables.  Then,  for  the 
random  vector  parameters  case,  the  CRLB  can  be  derived  by  a 


straightforward  modification  of 

the 

der ivat ion 

for 

the 

non-random  parameter  cases.  From 

the 

der ivat ion 

of 

the 

non-random  parameter  case,  we 

can 

easily  understand 

the 

fundamental  concept  of  the  CRLB 

and 

the  basic 

relation 

between  the  CRLB  and  the  variance  of  the  estimation  error. 

6.4.1  CRLB  for  Non-random  Variable  Case 

First,  we  assume  X is  an  unknown  constant  parameter  to 
be  estimated  from  a sequence  of  measurements  y^,y2»...»yK 

rp  a 

as  shown  in  Fig.  6.1  where  (Y  i»Y  2'  • * * »Y  k'  * Assuming  X 
is  an  unbiased  parameter  estimate  we  have 

Xf(X|X)dX  = X (6.12) 

or  from  Eq.  (6.1),  we  have 
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n 


X 


(6.13) 


II  ...j  T(y)f(y|X)dy  = 


K-fold 

integral 

A 

where  X=T  (^)  . Differentiating  both  sides  of  Eq.  (6.13) 
with  respect  to  X we  have 


rr  ( 3f (y |X) 
J---jT(y)  9 x-dy 


X ^ 


K-fold 

integral 

Now  we  rewrite  Eq.  (6.14)  as 


(6.14) 


r3f  (y|X)-| 


integral 


[•••  [T (y ) 

K-fold 

integral 


9£nf  (y  | X) 

) yy—  If  (y|X)dy  = 1 


Inspection  of  Eq.  (6.16)  shows  that  it  is 

9Jlnf  (y|X) 

E(T(y) ^ ] = 1. 


(6.15) 


(6.16) 


(6.17) 


Next,  we  examine  the  normalized  correlation  coefficient 

^ 9£.nf(yjx) 

between  X=T  (^)  and  §X  • From  the  definition  of P , we 
have 


3ilnf(y|X)  3 P,nf  (y  | X) 
E{  [ (T(y)-E(T(y)  ) ] [ ^ E( ) 


/var~fTy~)- Vvar(”n^^lX  ) 


where  var  denotes  variance.  Now, 

enf  (^|  X)  1 f 1 3 f (y  | X)  "I 

E [ Jx  | = E jFTyJx)  5x  J 
rr  r f<zix>  i 

= JJ-'-j  — sx-  • nzpo  f <nlx>di: 

- A jj  ■■■  |f<r lx>dn  - w ui  - o- 


6.18 


(6.19) 


hence 


E[4t  inf  (y | X)  ] = 0. 


(6.20) 


From  Eq.  (6.20)  and  simple  manipulations,  Eq.  (6.18) 


becomes 


E[T(y)^nf(y_X)_] 


/varT(y)^var(3iLnf  (yTxT7 


(6.21) 


From  Eq.  (6.17) , then 


/ st7 — r / 3Hnf  (y  X) 

'varT(il  -yvar( , 


(6.22) 


From  the  definition  of  p,  we  know  that  p is  equal  to  or 
less  than  1.  Therefore,  we  obtain 





var  (T(y)  ) = var  X >_ 


3 X 


(6.23) 


or  equivalently 


var  X = E [ (X-X) ] > Stnf[  |x)  2 . (6.24, 

" 3X 1 1 

Equation  (6.24)  is  called  the  Cramer-Rao  inequality  for  the 
unbiased  estimate  [6-15].  Note  that  the  CRLB  is  a bound  on 
the  mean-square  error.  For  a biased  estimate,  the  CRLB  is 


var  X > 


4 , dB (X) 
dX 

E((B"f '*!.*>  ,2)' 


(6.25) 


where  B(X)  is  the  bias  function  of  X. 


For  the  non-random  variable  vector  case,  we  have 
directly 


var  [x^-x^]  J 


(6.26a) 


where  J11  is  the  iith  element  in  the  KxK  square  matrix  J-1. 


The  elements  of  J are 


, a ,r8t.np(E 

ij  EL  SXj 


I 3Hnp(y.lx) 


lx)1 

j J 


X-  [x1,x2 , . . . ,xNl  , £=  [y1,y2, . . . ,yNl 


^6. 26b) 


where  x is  a vector  of  non-random  variables,  to  be 


estimated  and  ^ is  the 
Equ ivalently , 


measurement  data  vector. 


where 


J = E{  (Vx  Unp  (y  I x)  ] ) (Vx  Unp  (y  | x)  ] ) 1 } 


y k r_l_  JL_  _A_iT 

vx  9x  ' 9x_ ' ’ * ‘ ' 9x.  J 


(6.27) 


k 

The  J matrix  is  commonly  called  Fisher's  information  matrix 
16-1,  6-2,  6-3,  6-4,  6-15]. 


6.4.2  CRLB  for  Random  Variable  Vector  Case  [6-1 ] . 


I 


4 


For  the  random  variable  vector  case,  the  information 


matrix  JT  now  consists  of  two  parts 


JT  a JD  + Jp 


(6.28) 


where 


JD  = E({VxUnp(y|x)  ] }{VxUnp(y|x)  ] }T)  , (6.29) 


and 


Jp  = E ( { Vx^np (x) ) { Vxdnp (x) }T) . 


A „ f9£np  (y.[  x)  . 9Hnp(y,[x)  ] 

D.  . ~ b 9x.  9x . ' 

i]  L i J 

\ _["9!lnp(x)  9Hnp(x)"I 

ij = EL~^r  * -^rJ 


(6.30) 


il 

The  matrix  JD  represents  information  obtained  from  the  data 
or  from  the  probability  density  P(^|x)  of  the  MAP  estimate. 
The  matrix  Jp  represents  information  obtained  a pr ior i . 
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The  correlation  matrix  of  the  error  is 


< 

I 


i 


« 


n A „ , T. 

Re  = (6.31) 

A 

where  x£=  (x-x) . The  diagonal  elements  represent  the 
mean-square  errors  and  the  off  diagonal  elements  are  the 
cross  correlations.  The  mean-square  error  of  the  estimate 
as  a function  of  the  information  matrix  is 


Elxe  1 1 (J^1)11. 

l 

The  diagonal  elements  in  the  inverse 
information  matrix  JT  are  the  lower 
corresponding  mean-square  errors,  and  this 
of  interest  here. 


(6.32) 

of  the  total 
bounds  on  the 
is  the  situation 


6.5  Derivation  of  the  CRLB  for  MAP  Estimates  with  a 
Poisson  Noise  Model 

The  estimate  of  error  covariance  is,  in  general,  very 
complicated  to  find  due  to  the  complexity  of  the  posterior 
density.  However,  for  the  MAP  estimate  it  is  possible  to 
derive  an  expression  for  a lower  bound  on  the  variance 
because  we  know  the  a pr ior i density  p(^|x)  and  probability 
density  p(x)  of  x.  From  Eq.  (6.28),  we  have 


J 


T 


+ J 


P' 


(6.28) 


where  JD  and  Jp  are  defined  in  Eqs.  (6.29)  and  (6.30) 
respectively.  From  the  last  chapter,  we  have 
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l 


i 


NJO'  *.  V r- 


p (d | f ) = n 


d./a  -Xb. 
(Abi)  e 

cf 

a(^)! 


(6.33) 


p(f)  = Kbexp{-2-(f-f)TR^1(f-f)  } . 


(6.34) 


From  Eqs . (6.30)  and  (6.34),  we  can  obtain  [6-1,  6-10, 


6-11] 


JP  = RfX 


(6.35) 


where  Rf  is  the  covariance  matrix  of  the  image.  From 
Eqs.  ( 6 . 2!> ) and  (6.33),  we  have 

JD  = El  [Vf  inp  (d  | f ) 1 [VfHnp(d|  f)  1 ) 


= El (XHT(q-l) ] [XH  (q-1) ] M 
= \2HTE[  (q-1)  (q-l)TlH, 


(6.36) 


where  E denotes  expectation.  From  Eq.  (6.10)  we  have 


E[q]=l,  thus. 


where 


Jp  = X2HTE  [ (q-q)  (q-q)  T]  H 


2 T 

= ( H R H 

D a 


R^  = El (q-q) (q-q)T) 


(6.37) 


(6.38) 


Substituting  Eqs.  ^6.38)  and  (6.35)  into  Eq.  (6.28),  we 
have 


2 T - 1 

JT  = \ H R^H+Rf 1 


(6.39) 


This  is  the  total  information  matrix  JT  for  the  MAP 
estimate  with  the  Poisson  noise  model. 

When  the  J-p1  exists,  from  Eq.  (6.32)  we  have 


(6.40) 


where  (JT1)  are  the  diagonal  elements  of  j"1  and  f^  is 


an  estimate  of  the  ith  component  of  the  restored  image 
vector  f . Inspection  of  Eq.  (6.40)  indicates  that  the 
error  bound  depends  on  four  quantities:  the  Rate  function 
constant  X,  the  discrete  blurring  matrix  H,  the  covariance 
matrix  R^,  and  the  covariance  matrix  of  the  image  R^ . To 
obtain  some  physical  meaning  from  this  expression,  we 
assume  that 


R = o I 

a a 


and 


Rf  = ofI 


to  obtain 


.-1 


( (X2HTHo^I+(o^)"1I) )_1 
= I (a2l|H||2o2I)  + (o2)"1I]‘1  , 


(6.41) 


Thus 


A 9 f 


+ 1 


(6.42) 


From  Eq.  (6.42),  we  divide  out  (ij  on  both  the  numerator  and 
the  denominator  of  the  right  hand  side  of  Eq.  (6.42)  and 
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z 


* 


get 


E[(f  -f  )2]  = -j Jr-  - 

1 1 SfH 

°f 

From  Eq.  (6.43),  we  observe  that 


(6.43) 


2 2 

(1)  When  a ^ is  large,  we  can  drop  the  term  l/of  from  the 


denominator , 


Multiplying  both  the  numerator  and 

2 


denominator  of  Eq.  (6.43)  by  the  mean  intensity  5^,  gives 


2 

E [ ( f i~f i ) 1 = 


A 


:xbi)2|(H||2o2 


(6.44) 


In  this  situation,  , H,  and  play  more  important  roles 

in  the  error  bound  of  the  Poisson  noise  model  than  the 

2 

variance  of  the  object  a_; 

(2)  The  error  bound  decreases  with  the  square  of  the 
ensemble  mean  rate  function  (Xb^).  From  the  experimental 
results  of  Fig.  6.2,  we  see  that  increasing  the  ensemble 
mean  rate  function  of  the  Poisson  process  Xb^,  reduces  the 


Poisson  noise  degradation.  Consequently,  when  the  (SNR)rms 
is  greater  than  or  equal  to  10  db  (Xb^lOO),  then 


E(  (fi-fi)2)  > 10"4  , 


(6.45) 


and  the  Poisson  noise  degradation  effects  are  small  in  a 
practial  sense. 

(3)  The  error  bound  is  also  inversely  proportional  to  the 
squared  norm  of  the  point  spread  function  H and  the  error 
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Images  with  Poisson  noise  at  different 
(SNR»rms 

(A)  Original  object  image 

(B)  The  Poisson  noisy  image  with 
(SNR)rms./5 

(C)  The  Poisson  noisy  image  with 
(SNR)rms=/I° 

(D)  The  Poisson  noisy  image  with 


1 


due  to  noise  is  "amplified"  by  the  point  spread  function. 
This  results  from  the  ill-conditioning  of  the  restoration 
process. 

6.6  Conclusions 

We  conclude  that  the  MAP  estimate  for  the  Poisson 
noise  model  is  an  unbiased  estimate  and  have  found  the 
Cramer-Rao  lower  bound  (CRLB)  for  the  variance  of  the 
estimation  error.  The  CRLB  is  the  tightest  lower  bound  for 
an  efficient  estimate.  When  an  efficient  estimate  does  not 
exist,  the  lower  bound  can  be  improved  compared  to  the 
Cramer-Rao  inequality.  Better  lower  bounds  may  be  the 


Bhattacharyya  bound 

and 

Barankin 

bound,  but  these 

bounds 

are  very  difficult 

and 

tedious 

computationally. 

The 

Barankin  bound  does  not  require  the  probablity  density  to 
be  differentiable  and  it  gives  the  greatest  lower  bound. 

It  requires  a maximization  over  a function  to  obtain  the 
bound  and  the  procedure  for  finding  this  maximum  is  usually 
not  straightforward. 

Use  of  the  statistical  estimation  method  is 
particularly  desirable  from  the  viewpoint  of  error  analysis 
because  known  techniques  can  be  applied  to  compute  the 
error  bound.  We  have  developed  the  CRLB  for  MAP  estimation 
with  the  Poisson  noise  model  and  also  shown  the  behavior  of 
the  CRLB  approximation.  From  these  facts,  we  are  able  to 
find  the  algoritm  whose  mean-squared  errors  is  closest  to 
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the  CRLB 


It  must  be  remembered  that  the  CRLB  is  a lower 


bound  and  that  the  actual  restoration  error  will  be 
greater.  It  is  possible  that  a better  suboptimal  sectioned 
MAP  algorithm  can  be  found  to  reduce  the  actual  restoration 
error  closer  to  the  Cramer-Rao  lower  bound. 


CHAPTER  7 


COMPARISON  BETWEEN  THE  LMMSE  RESTORATION  FILTER 
AND  THE  MAP  RESTORATION  FILTER 

7. 1 Introduction 

The  most  common  goal  for  restoring  a degraded  image  is 
to  reduce  the  measurement  error  and  to  make  the  information 
more  visible.  A true  comparison  between  filters  should 
follow  some  objective  criteria.  However,  the  mechanism  of 
human  information  extraction  is  not  well  understood  and 
there  are  no  universally  agreed  upon  criteria  by  which  to 
judge  the  quality  of  a proposed  image  restoration  filter. 
There  are  two  simple  criteria  which  are  commonly 
accomodated.  One  is  the  numerical  closeness  of  the 
restored  image  to  the  undegraded  original  object  image  in 
terms  of  mean  square  error,  and  the  other  is  the  visual 
subjective  appearance  of  the  restored  image  compared  to  the 
original.  These  two  criteria  are  often  in  conflict  because 
of  the  complex,  nonlinear,  and  adaptive  properties  in  the 
psycho-physical  processes  of  human  vision  [7-1].  Pearlman 
has  proposed  a new  compromise  criterion  [7-2]  which  is  a 
weighted  mean  square  error.  The  amplitude  weights  are  not 
constants  but  are  dependent  upon  the  contrast  ratio  of  the 
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image  and  also  upon  an  exponential  function  of  average 


intensity.  This  new  weighted  mean  square  error  criterion 
reflects  some  degree  of  the  nonlinearity  and  complex 
properties  of  the  human  visual  system,  although  it  is 
complicated  to  evaluate.  The  numerical  closeness  criterion 
is  most  often  employed  because  it  is  well  defined  and 
mathematically  tractable.  There  is  a fundamental 
difference  in  the  estimation  criterion  between  the  MAP  and 
the  LMMSE  image  restoration  filters.  The  LMMSE  filter  is 
based  upon  minimizing  linear  minimum  mean-square  error 
under  perfect  a priori  knowledge  of  the  object  and  image, 
while  the  MAP  filter  is  based  upon  the  maximization  of  the 
a poster ior i density  of  the  image.  Although  it  is 
theoreticaly  hard  to  make  any  comparison  between  them 
because  of  this  fundamental  criterion  difference,  we  will 
compare  them  based  upon  the  first  two  quality  criteria  and 
computation  time. 


The  structure  of  the  LMMSE  filter  is  discussed  and 
analyzed  in  the  next  section.  In  section  7.3  we  discuss 
the  implementation  and  illustrate  additional  experimental 
results.  In  section  7.4  we  discuss  image  quality  measures 
and  compute  the  normalized  mean  square  error  of  the 
restorations.  In  the  final  section  we  make  comparisons  and 
state  the  conclusions. 


> . ;c*j4n 


7.2  The  Structure  of  the  LMMSE  Filter  and  MAP  Filter 


7.2.1  Structure  of  the  LMMSE  Restoration  Filter 

From  chapter  3,  the  LMMSE  filter  ^transfer  function 
W(u,v)  is  derived  by  minimizing  e(  //•2  (x,y)dxdy).  By 

— OO 

Parseval's  theorem,  it  is  equivalent  to  minimizing 

oo 

E [JJ  | e(u,v)|  audv]  , where  e is  the  Fourier  transform  of 

— OO 

e(x,y).  Use  of  this  minimizing  mean-square  error  criterion 
and  the  orthogonality  principle  [7-10,  7-11,  7-12]  yields 
the  transfer  function  of  the  LMMSE  filter  as 


Wp(u,v) 


gg 


( u , v ) 
(u,v) 


(7.1) 


where  4^  is  the  cross-spectral  density  of  the  detected 
image  g(x,y)  and  object  f(x,y)  while  ‘J’gg  is  the  spectral 
density  of  the  detected  image  g(x,y).  From  a straight 
forward  substitution  into  Eq.  (7.1)  [7-3, 7-4],  we  have 


NV*  (u,v)  4>f  (u,v) 

W ( u , v ) = 5 

p l+N||V(u,v)|r4>f  (u,v) 


(7.2) 


where  N is  mean  number  of  photon  counts  in  the  detected 
image  V(u,v)  is  the  Fourier  transform  of  the  PSF  h(x,y)  and 
$f(u,v)  is  the  spectral  density  of  the  object.  The 
detailed  derivation  of  Eq.  (7.2)  is  described  in  [7-3]. 
For  a linear  Gaussian  additive  noise  model,  the  Wiener 
filter  is  [7 -9,  7-13,  7-14,  7-15,  7-16] 
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Ww(u,v)  = 


V*  (u,  v)  <t>f  (u,  v) 


V(u,v)||  4>f  (u,  v)  +$n  (u,v) 


(7.3) 


where  <t>  n(u,v)  is  the  spectral  density  of  noise  which  is 
statistically  independent  of  the  signal,  ^(UjV)  is  the 
spectral  density  of  the  object  and  ^(u,v)  is  the  Fourier 
transform  of  h(x,y).  Rewriting  Eq.  (7.3),  the  Wiener 
filter  takes  the  most  familiar  form 


where 


„ , 4 V*(u,v) 

W (u,v)  = '-j—y 

W ||«(u,v)||2+| 


Vu,v) 


(7.4) 


S is  called  the  signal-to-noise  (SNR)  for  the  linear 
additive  Gaussian  noise  model. 


Since  the  linear  additive  Gaussian  model  assumes 
signal-independent  noise,  it  is  reasonable  to  define  B as 
the  signal-to-noise  ratio.  Similarly,  from  Eq.  (7.2)  for 
the  Poisson  noise  model,  the  LMMSE  filter  is 


& * / II  V ) 

W ( u , v ) = -^4  r 

P ||  V (u,v)||2+^ 


(7.5) 


where  a=N  $f(u,v)  . The  function  Wp(u,v)  is  the  same  as 
W^u,v)  except  that  a is  defined  differently  from  B.  since 
Poisson  noise  is  signal-dependent,  the  SNR  is  not  well 
defined.  However  a can  be  called  the  equivalent 


n 


■ :>■  • - Uti i'i  ..  - 7*J 


signal-to-noise  ratio  (SNR)  . From  Eq.  (7.5),  when  the 

eq 

value  of  a is  very  large,  which  is  the  case  when  the  rate 
function  is  high,  then  Eq.  (7.5)  becomes 


W (u,v)  ~jy-1(u,v) 


(7.6) 


Thus,  the  LMMSE  filter  approaches  the  inverse  filter 
in  the  absence  of  Poisscn  noise.  Indeed,  the  larger  the 
rate  function,  the  lesser  is  the  degradation  due  to  Poisson 
noise.  In  this  case,  the  LMMSE  filter  only  needs  to  remove 
the  blurring  degradation  effects.  As  discussed  in 
chapter  5,  Poisson  noise  effects  are  much  more  pronounced 
at  low  light  levels  when  the  value  of  a is  smaller.  In 
this  case,  the  LMMSE  filter  Wp  is  dominated  by  Poisson 
noise  and  image  signals  will  be  seriously  distorted.  Thus 
the  performance  of  the  LMMSE  filter  will  expected  to  be 
worse  at  lower  equivalent  SNR's.  Although  the  LMMSE  filter 
is  based  upon  the  minimum  mean-square  error  criterion,  this 
estimation  error  is  a minimum  under  the  condition  that  the 
a pr ior i knowledge  is  perfect.  Functions  such  as  the 
spectral  density  of  the  object  and  the  mean  number  of 
photon  counters  N must  be  perfectly  known.  In  reality,  $£ 
and  N are  never  perfectly  known  and  must  be  estimated  from 
the  observation  [7-7,  7-18].  Hence,  the  actual  LMMSE  error 
does  not  reach  the  minimum.  In  short,  the  LMMSE  filter 
tries  to  force  the  solution  toward  the  inverse  solution 

with  some  sort  of  smoothness  controlled  by  the  equivalent 
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7.2.2  Structure  of  the  MAP  Restoration  Filter 

It  has  been  noted  in  preceding  chapters  that  the 
fundamental  MAP  estimate  contains  maximum  likelihood  (ML) 
and  a priori  terms  in  its  solution.  From  previous 
chapters,  the  MAP  estimate  equation  is 


AHT(q-l)-R^1(f-f)  = O , 


(7.7) 


or  equivalenty 


Imap  = i+XRfH  <s-P 


(7.8) 


The  physical  interpretation  of  the  MAP  estimate  is  that 

maximizing  the  probability  p ( d | jf ) forces  the  solution 

toward  the  inverse  solution  which  is  the  maximum  likelihood 

solution,  while  maximizing  the  probability  p (JE)  is 

equivalent  to  enforcing  a smoothness  criterion.  Thus,  the 

MAP  filter  tries  to  balance  the  inverse  solution  with  a 

smoothness  constraint  (7-5] . Another  physical 

interpretation  from  Eq.  (7.8)  is  that  the  MAP  estimate 

tries  to  move  the  solution  of  the  estimate  f_._  from  the  a 

~ —MAP  — 

pr  ior  i nonstationary  mean  jE  to  a maximum  likelihood 

solut ion  f T . 

—ML 

7.2.3  Conclusions 


The  LMMSE  filter  and  MAP  filter  are  performing  similar 


1 


4 
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functions 


of  balancing  the  inverse  solution  with  a 


smoothness  constraint. . The  LMMSE  filter  uses  the 
equivalent  SNR  to  control  the  balancing,  while  the  MAP 
filter  uses  the  covariance  matrix  R of  the  object  as  a 
measure  of  the  confidence  in  the  nonstationary  mean  jf  and 

A 

the  maximum  likelihood  solution  jfML  as  a solution  to  the 
restoration  filter.  The  LMMSE  filter  is  based  on  the 
assumption  that  the  detected  image  intensity  can  be 
approximately  modeled  by  a stationary  random  field.  For  a 
typical  image,  each  part  of  an  image  generally  differs 
sufficiently  from  other  parts  so  that  the  stationarity  is 
not  generally  valid.  Moreover,  the  LMMSE  filtering  process 
is  insensitive  to  abrupt  changes  of  image  signals.  This 
results  in  edge  smoothing  and  contrast  reduction.  The  MAP 
filter  does  consider  the  nonstationary  properties  of  random 
image  fields  because  it  contains  an  ML  term  and  an  a pr ior i 
term.  Also  it  is  an  adaptive  processor  which  depends  on 
the  local  nonstationary  properties  of  image  signals.  The 
MAP  filter  theoretically  needs  more  a pr ior i knowledge,  but 
in  actual  implementation,  the  MAP  filter  uses  less  a pr ior i 
knowledge  than  the  LMMSE  filter.  Both  filters  require 
knowledge  of  the  blurring  matrix  H and  mean  number  of 
photon  counts  fT  but  the  LMMSE  filter  also  requires  the 
spectral  density  of  the  object  It  can  be  intuitively 
concluded  that  the  MAP  filter  should  perform  better  than 
the  LMMSE  filter  for  the  Poisson  noise  model. 
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7.3  Implementation  and  Experimental  Results  of  the  LMMSE 
Filter  and  the  MAP  Filter 


Our  experimental  implementation  of  the  LMMSE  filter  is 
based  on  Eq.  (7.5)  using  a fast  Fourier  transform 
algorithm.  The  ensemble  mean  photon  counts  N and  the 
object  spectral  density  4> ^ are  estimated  from  observed 
photon  count  data.  The  estimate  of  can  be  made  by 
substituting  a similar  "prototype"  spectral  density 
suggested  by  Cannon  in  the  blind  deconvolution  process 
17-6]  , or  estimated  by  an  iteration  method  suggested  by 
Limb  in  a method  of  image  restoration  called  spectral 
subtraction  (SS1R)  17-7].  Unfortunately,  they  assume  the 
noise  is  linear  signal-independent  additive.  Because  the 
Poisson  noise  is  signal-dependent,  the  estimation  of  is 
different. 


The  spectral  density  of  the  object  ^ is  related  to 
the  spectral  density  of  the  detected  image  by 


<J>d(u,v)  = N+  (N)  2 ||v(u,v)||  2<)>f  (u,v) 


(7.9) 


where  N is  ensemble  mean  number  of  photon  counts  and  V(u,v) 
is  the  Fourier  transform  of  the  PSF  h(x,y).  Thus,  the 
spectral  density  of  the  object  ^ can  be  estimated  by  an 
iterative  method,  although  this  involves  an  inverse  filter 
with  V.  we  do  not  investigate  the  estimation  of  the 
spectral  density  $£.  Instead  we  assume  that  is  a white 
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spectral  density  in  implementat i ng  the  LMMSE  filter.  The 
white  object  spectral  density  assumption  extracts  more 
higher  spatial  frequency  information  such  as  edge  and  fine 
uetaii  of  the  image  because  the  spatial  frequency  content 
of  many  images  falls  rapidly  at  high  spatial  frequencies. 


For  generality,  a two-dimensional  moving  average 
blurring  point  spread  function  (PSF)  is  chosen  for  the 
simulation  rather  than  a Gaussian  blurring  PSF  because  it 
has  singularities  and  phase  reversals  in  the  frequency 
response.  Restored  images  with  the  LMMSE  filter  for 
diferent  equivalent  SNR's  are  illustrated  in  Fig.  7.1. 
The  amount  of  restoration  is  controlled  by  the  a.  The 
serious  effects  of  ill-conditioning  can  be  seen  in 
Fig.  7.1.  The  restored  image  gradually  blows  up  as  a goes 
to  higher  values,  while  it  becomes  more  noisy  as  a goes  to 
lower  values.  Because  the  Poisson  noise  degradation  is 
very  pronounced  for  small  values  of  a and  the 
i 1 1-cona i t ioning  of  the  inverse  filter  is  worse  at  larger 
values  of  rt,  implementation  of  the  LMMSE  filter  is  very 
sensitive  to  the  equivalent  SNR  a. 


Figures  7.2  and  7.3  show  results  with  the  nonlinear 
sectioned  MAP  filter  implemented  with  Newton-Raphson 
iterative  techniques  as  before.  Figure  7.2  has  results 
with  no  blur,  and  Figure  7.3  shows  results  with  linear 
blur.  The  experiments  for  both  cases  of  degradation  are 
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Figure  7.1  The  restored  images  with  the  LMMSE 

filter  for  different  estimated  (SNR)  6 

eg 
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Images  restored  by  the  LMMSE  filter  and 
the  MAP  filter  at  (SNR)  =/2 . 5 


Original  object  image 
Poisson  noisy  image 
Image  restored  by  the  LMMSE  filter 
Image  restored  by  the  MAP  filter 


Figure  7.2b  Images  restored  by  the  LMMSE  filter  and 

the  MAP  filter  at  (SNR)  = /5 

rms 


(A)  Original  object  image 

(B)  Poisson  noisy  image 

(C)  Restored  image  by  the  LMMSE  filter 

(D)  Restored  image  by  the  MAP  filter 
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Figure  7.2c  Images  restored  by  the  LMMSE  filter  and 

the  MAP  filter  at  (SNR)  =/T0 

rms 

(A)  Original  object  image 

(B)  Poisson  noisy  image 

(C)  Restored  image  by  the  LMMSE  filter 

(D)  Restored  image  by  the  MAP  filter 


7. 2d  Images  restored  by  the  LMMSE  filter  and 
the  MAP  filter  at  (SNR)  =/20 


Original  object  image 
Poisson  noisy  image 

Restored  image  by  the  LMMSE  filter 
Rer“-ored  image  by  the  MAP  filter 


Images  restored  by  tne  uim&c.  rater  arm 
MAP  filter  with  two-dimensional  linear 
blurring  degradation  at  (SNR)  rms='2" 5 

(A)  Original  object  image 

(B)  Degraded  image 

(C)  Restored  image  by  the  LMMSE  filter 

(D)  Restored  image  by  the  MAP  filter 


Images  restored  by  the  LMMSE  filter  and  the 

MAP  filter  with  two-dimensional  linear 

blurring  degradation  at  (SNR)  = /5 

rms 

(A)  Original  object  image 

(B)  Degraded  image 

(C)  Restored  image  by  the  LMMSE  filter 

(D)  Restored  image  by  the  MAP  filter 


Images  restored  by  the  LMMSE  filter  and  the 
MAP  filter  with  two-dimensional  linear 
blurring  degradation  at  (SNR)  rms=/TfT 

(A)  Original  object  image 

(B)  Degraded  image 

(C)  Restored  image  by  the  LMMSE  filter 

(D)  Restored  image  by  the  MAP  filter 


Figure  7.3d  Images  restored  by  the  LMMSE  filter  and  the 
MAP  filter  with  two-dimensional  linear 
blurring  degradation  at  (SNR)  rms=/2"0 

(A)  Original  object  image 

(B)  Degraded  image 

(C)  Restored  image  by  the  LMMSE  filter 

(D)  Restored  image  by  the  MAP  filter 
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performed  with  different  (SNR)rms. 

7 . 4 Numer ical  Values  of  Restored  Image  Quality  Measures 


One  of  the  most  commonly  used  numerical  quality 
measures  is  the  mean  square  error  (MSE) . It  can  be  defined 
as 


N N 

MSE  = £ £ [f (m,n) -f (m,n) ] 
m=l  n=l 


(7.10) 


Although  this  measure  is  attractive  because  it  is 
tractable,  it  does  not  match  human  evaluation  on  many  types 
of  images.  It  is  also  possible  to  define  a measure  based 
on  the  MSE  and  energy  normalization  (7-9).  This  normalized 
mean  square  error  (NMSE)  is  defined  as 


NMSE  = 


N N 

^ £ [f (m,n)-f (m,n) ] ' 

m=  1 n=l 

N N , 

£ £[f(m,n)P 

m=l  n=l 


(7.11) 


The  NMSE  performs  somewhat  better  than  MSE  and  retains  the 
analytic  tractability.  For  these  reasons,  we  use  it  in 
this  chapter.  Other  numerical  measures  such  as  normalized 
error  (NE)  , and  Laplacian  mean  square  error  (LMSE) 
[7-af7-9]  etc.  will  not  be  used  here.  Based  on 
Eo.  (7.11),  the  NMSE  of  the  restored  image  between  the 
LMMSE  filter  and  the  MAP  filter  will  be  computed.  The 
results  are  shown  in  Table  7.1  and  Table  7.2  for 
non-blurring  and  blurring  cases  respectively. 
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7 . 5 Compar isons  and  Conclusions 

The  image  quality  of  the  pictures  restored  by  the  MAP 
filter  shown  in  Figs.  7.2  and  7.3  seems  superior  to  those 
processed  by  the  LMMSE  filter.  Images  restored  by  the 

LMMSE  filter  have  an  excessive  enhancement  of  Poisson 

/\ 

noise,  especially  for  higher  values  o fa.  in  addition,  the 
NMSE  of  the  MAP  filter  is  lower  than  that  of  the  LMMSE 
filter.  The  LMMSE  filter  strongly  depends  on  perfect  a 
pr ior i knowledge  of  the  object  and  it  is  very  sensitive  to 
parameters  in  the  filter  implementation.  Table  7.1  and 
Table  7.2  show  that  the  MAP  filter  has  advantages  over  the 
LMMSE  filter,  particularly  at  lower  (SNR)rms . Furthermore, 
the  NMSE  of  the  LMMSE  filter  does  not  improve  greatly  at 
increasing  (SNR)rms  as  it  does  with  the  MAP  filter  in  the 
blurring  case.  A possible  explanation  for  this  is  that  the 
LMMSE  filter  performs  very  well  for  the  higher  (SNR),  but 
the  MAP  filter  also  works  better  for  the  lower  (SNR)  at 
which  Poisson  noise  dominates.  The  disadvantage  is  that 
the  MAP  filter  is  a nonlinear  spatial  estimate  which  needs 
iterative  methods  for  solution.  Thus,  the  computing  time 
of  the  MAP  filter  is  longer  than  that  of  the  LMMSE  filter. 
The  better  performance  is  achieved  at  the  expense  of 
additional  computing  time. 
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Table  7.1 


The  NMSE  for  nonblurrina  cases 


(SNR) 

rms 

— 

LMMSE  filter 

MAP  filter 

f2~.  5 

0. 3601287 

0.15927 

/5 

0.18131 

0.11352 

/To 

0.090158 

0.072803 

/2  0 

0.0449264 

0.04180208 

Table  7 . 2 

The  NMSE  for  blurrinq  cases 

(SNR) 

rms 

LMMSE  filter 

MAP  filter 

/2T5 

0.269397 

0.190431 

/5 

0.1502184 

0. 0994055 

/To 

0.1466819 

0.054126 

/20 

0.1433148 

0.036138 
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SUMMARY  AND  CONCLUSIONS 


In  this  work  we  have  modeled  photon  resolved  image 
signals  as  a Poisson  point  random  process  and  developed  an 
optimal  spatial  restoration  filter  for  the  Poisson  noise 
model . 

Poisson  noise  is  an  inherent  part  of  any  detected 
image  and  is  particularly  evident  in  low  level  image 
signals.  Because  it  results  from  the  discrete  random 
nature  of  quantum  limitations,  it  is  signal-dependent.  The 
optimal  spatial  filter  was  based  on  a criterion  of 
maximizing  the  a poster ior i probability  density.  The 
formulation  and  solution  of  the  MAP  estimation  problem  have 
been  presented.  It  has  been  found  that  the  overlap-save 
sectioning  method  with  a Newton-Raphson  iterative  solution 
is  the  most  efficient  way  of  coping  with  the  nonlinearity 
and  large  dimensionality  of  the  MAP  estimation  equations. 
The  implementation  of  the  MAP  filter  with  the  Poisson  noise 
model  was  made  for  both  blurring  and  non-blurring 
degradation  cases.  It  has  been  demonstrated  that  the  MAP 
filter  with  the  Poisson  noise  model  has  improved 
performance  because  the  MAP  filter  can  be  generalized  to 
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linear  or  nonlinear  image  models  and  to  noise  models 
different  from  additive  Gaussian  noise.  In  addition,  the 
MAP  filter  can  be  a local  adaptive  processing  filter  and 
can  be  extended  to  space-variant  blurring.  It  also  has 
been  shown  that  modeling  images  with  a nonstationary  mean 
and  stationary  variance  gives  useful  a priori  information 
for  the  MAP  filter. 

The  Cramer-Rao  lower  bound  (CRLB)  on  the  mean-square 
estimation  error  of  the  MAP  unbiased  estimate  was  derived 
for  the  Poisson  noise  moael.  It  is  likely  that  the  CRLB 
may  be  useful  for  finding  the  best  suboptimal  sectioned  MAP 
filter.  A comparison  was  made  between  the  LMMSE  filter  and 
the  MAP  filter  with  the  Poisson  noise  model.  It  has  been 
shown  that  the  quality  of  the  restored  image  of  the  MAP 
filter  is  superior  to  that  of  the  LMMSE  filter  by  simple 
subjective  evaluation  and  by  numerical  closeness  criteria. 

Boulter  (b-lj  has  shown  that  even  with  large  amounts 
of  blurring  present  in  a detected  image,  a low  noise  level 
permits  almost  complete  restoration.  However,  photon 
resolved  image  signals  have  a very  low  s ignal-to-noise 
ratio,  making  it  impossible  to  obtain  perfect  restoration 
for  image  signals  suffering  from  both  blurring  and  Poisson 
noise . 


The  research  pursued  in  this  dissertation  may  be 
extended  in  several  areas.  A more  detailed  study  of  an 
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optimal  local  adaptive  MAP  filter  using  the  local 
properties  of  the  first  and  second  moment  of  image  signals 
would  be  of  considerable  interest.  Use  of  these  local 
statistics  is  interesting  because  image  fields  are 
inherently  nonstationary.  Another  area  of  practical 
importance  is  to  develop  a fast  algorithm  for  the  MAP 
filter  or  a recursive  MAP  filter  for  saving  computing  time 
and  memory  space.  A recursive  MAP  filter  would  not  only 
offer  computational  advantages  over  a non-recursive  filter 
but  also  could  be  applied  to  space-variant  and 
nonstationary  models.  However  it  is  expected  that  the 
recursive  spatial  restoration  filter  would  be  very 
sensitive  to  errors  in  the  knowledge  of  the  PSF.  Another 
intereting  area  is  to  apply  these  results  to  other  types  of 
photon  resolved  image  signals  such  as  nuclear  medicine, 
medical  images,  astronomical  images  and  the  projection 
reconstruction  image  signals  (8-3].  In  these  cases,  the 
image  formation  system  model  and  data  acquisition  system 
must  be  carefully  modeled  to  identify  the  parameters  of  the 
photon  counting  system.  Another  interesting  possibility 
for  the  future  is  to  use  Lebedev's  "composite"  image  model 
to  develop  a "multicategory"  spatial  MAP  filter  (8-4,  8-5]. 
This  model  involves  multiple  categories  of  random  fields  in 
the  image  with  each  category  distinguished  by  its 
covariance.  A Gaussian  probability  density  is  associated 
with  the  occurence  of  each  category  in  the  data.  With  this 
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model,  the  optimal  estimate  is  obtained  by  filtering  the 
measurement  data  for  each  category,  and  forming  the 
estimate  as  the  weighted  sum  of  all  the  filter  outputs. 
The  weights  are  the  a poster  ior i probabilites  that  the 
point  is  a member  of  the  respective  category.  The 
composite  image  model  locally  decomposes  the  image  signals, 
and  it  should  be  able  to  model  the  local  nonhomogeneous 
information  in  images. 
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POISSON  RANDOM  NOISE  GENERATOR  [2-25] 


For  the  simulation  of  photon  count  observation  data, 
we  need  an  algorithm  for  generating  a sequence  of  random 
numbers  from  a population  conforming  to  the  Poisson 
distribution  with  mean.  The  solution  is  to  make  use  of  a 
random  number  generator  which  returns  a random  variable  z 
having  the  uniform  distribution  h(z)=l  for  0£z£l  and  h(z)=0 
for  other  values  of  z.  If  we  form  x k=z Q, z z 2, • • • , z k as 
the  product  of  a sequence  of  k+1  such  random  variables, 
then  the  lowest  value  of  k which  first  cause  *k  to  be  less 
than  or  equal  to  e \ will  be  a random  variable  which  has  a 
Poisson  distribution  with  mean  X.  a proof  of  this  property 
follows.  If  y and  z are  independent  random  variables  with 
pdf  g (y)  and  h(z)  respectively,  it  can  be  shown  that  [2-17] 
x=yz  has  pdf 

f(x)  = j g(y)h(^)^-  . (A.  1 ) 


Let  fk(xk)  denote  the  pdf  of  xk.  Since  f ( xQ)  =h  ( x0)  =1  we 
have,  for  x^=XqZ 
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fl<xl)  * 


dxn 

x 9(x0)h(3T)l r1 
X1  u x0  0 


(A. 2) 


fl(V  = 


1*1* = -in  x. 

X1  x0  1 


(A. 3) 


By  induction,  it  can  be  shown  that  xk=xk_izk  has  pdf 


fk(xk}  " nrr  in  \ 


(A. 4) 


By  use  of  the  well  known  recursive  relation 


where 


Ik  = x in  x-k  Ik_^, 


in  x dx. 


the  probability  that  x.  is  less  than  e is 


Pr(xk<e-X)  = J ^ £k(xk)dxk 


Thus 


= -(~ky--  [xinkx-k  XR_ x ] 6 


Pr<*kie"X>  = TFT  +Pr(xk-lle'X) 


(A. 5) 


(A. 6) 


(A. 7) 


Using  the  same  method  to  calculate  Pr^xk-Lie  ) • we  then 
have 


„ , - Av  -AfAk  , Xk_1 

Pr(xk<e  ) = e + -j^TT 


+ + A + 1]  (A. 8) 
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ano  we  conclude  that  if  x^.le  but  x^_^>e  , then  k obeys 
the  Poisson  distribution  with  mean  X.  This  algorithm  turns 
out  to  be  a very  accurate  end  fast  way  to  generate  the 
photon  counts  with  given  mean  image  intensities.  This 
algorithm  is  also  available  in  the  IMSL  subroutine  package 
(GGPOSH ) 12-26]. 


APPENDIX  B 

NEWTON-RAPHSON  ITERATIVE  SOLUTION  METHOD  FOR  NONLINEAR 
MAP  ESTIMATION  EQUATIONS 


A detailed  derivation  and  procedure  for  the 
Newton-Raphson  iterative  method  is  described  in  [ 4— 4 ] . 
Here,  we  summarize  the  Newton-Raphson  solution  of  the  MAP 
estimate  equations.  From  Eq.  (4.14),  we  get  N MAP  estimate 
equations,  that  is 


gl(^)  = (fi“1) S- (fl"?l)+^(f2_f2)  " °'  (B.  la) 

1 (1+p2)  1 1 A 2 2 


g.(f)  = (?f-i)+^(fi.1-Ii_1)-J(fi-Ii)^(fi+1-ri+1)(B.: 


= 0, 


i = 2,3,4,... , N- 1 


9Nlf>  - - O-(B.lc) 

£N  A w j.  in  j.  (i  + pz)  n « 


where  B = — tV-  , r = 1+Pa  • -L-  and 

1 + P 1-p**  Of 

P is  the  correlation  coefficients  between  pixels 


f is  the  object  vector  to  be  estimated 
f is  the  nonstationary  mean  vector 


Once  the  problem  has  been  set  up  as  in  Eqs.(B.la)  to 
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(B.lc),  the  solution  procedure  is: 

Step  1:  Choose  initial  guess  of  object  vector 


I 


f(k)_f(0)  =[f  f flT 

- L irl ' r2 ' * ' • ' N* 


(B. 2a) 


or 


f (k)  _ f (0)  _ r J J , ,T 

f f (d1*d2* • • • *dNl 


(B. 2b) 


where  f is  the  nonstationary  mean  and  d is  the  noisy 
observation  data  vector.  The  choice  of  initial  values 
affects  only  the  convergence  rate  because  the  equations 
finally  converge  to  a unique  set  of  roots.  Here,  the 
superscript  k denotes  the  kth  iterative  step 


Step  2:  Solve  the  linear  system  Eq.  (B.3)  to  find  the 


(k) 


solution  vector  6(f  ) 


4>(  f 


(k)  (k)  = 


(k) 


(B.3) 


where 


. / f (k ) . _ 3gi,f(k) 

*ij ' Tf7(-  } ' 

J 3 


_ r ,f(k).  „ ,f(kK  _ (k).  ,T 


and  incremental  vector 


5(f(k))  = [61(f(k),62(f(k)),...,6N(f(k))lT 


From  Eq.  (B.l)  and  4> . . (_f^  ) (f  ^ ) » we  obtain 


diagonal  entries  are  larger  than  that  of  off  diagonal 
entries.  Thus,  Eq.  (B.3)  equations  will  be  converged. 

Step  2 for  solving  the  linear  Eq.  (B.3)  can  be 

summarizea  as  follows: 

(k) 

(a)  Substitute  i into  Eq.  (B.l)  and  (B.4)  to  find  the 

, t < * (kh  =,  a 

values  of  £(.£  ) and  £(£  ) 

(b)  Solve  Eq.  (B.3)  using  a linear  equation  subroutine  to 

(k) 

find  the  solution  _6(jf  ) 

Step  3:  Update  the  approximation  to  the  root  for  the  next 
iteration 


f (k+1)  = f<k)+6(f(k))  ( B. 6 ) 

Step  4:  Check  for  possible  convergence  to  a real  root  _f  by 
applying  the  test 

6i(f(k))  < e for  all  i.  (B*7> 

Step  5:  It  Eq.  (B.7)  is  true  for  all  i,  then  f^  ^ ^ is 
taken  to  be  the  root.  If  Eq.  (B.7)  fails  for  any  i,  then 
the  process  is  repeated  starting  from  step  2. 

In  short,  this  iterative  algorithm  converges  very 
fast,  usually  in  about  two  to  three  steps  in  our 


simulation. 
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